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ABSTRACT 


The process of initiation and quasi-static propagation 
of crack under generalized plane stress and plane strain 
conditions in an isotropic, linear elastic niaterial is studied. 
The strain energy density failure criterion is used for pre- 
dicting crack growth characteristics in which stable crack 
growth is simulated by predicting a series of crack growth 
steps corresponding to piecewise loading increments when 
material ahead of the crack tip absorb a critical amount of 
elastic strain energy density, a material constant. Unstable 
crack growth which may lead to global instability, is predi- 
cted when the strain energy density factor assumes a critical 
value, which is again a material constant. The whole crack 
growth process from initiation to unstable growth is thus char- 
acterized by two material constants. 

The stress analysis is performed by a combined boundary 
element method (BEM) and finite element method (FEM) to exploit 
the advantages of both methods. The coupling of multi-domain 
BEM and FEM is achieved by minimisation of an energy functional 
leading to symmetric boundary element 'stiffness' matrix so 
that it could be easily incorporated into an existing finite 
element software. Starting with an initial structural dis- 
cretization of the cracked body, the model automatically 
computes stress-intensity factors, direction of propagation 
and amount of propagation for a given loading condition. The 
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process is repeated, again automatically, until the crack 
instability occurs or required number of crack steps is reached 
'Traction singular' guarter-point boundary elements which 
contain the dominant l//r traction singularity and dominant /r 
displacement variation around the crack tip are used at the 
crack tips. 

The quasi-static crack propagation analysis system 
requires an 'interactive simulation' environment which 
makes effective use of recent advances in computers, display 
devices, automatic mesh generation, computer graphics and 
data base management systems. The finite element-bomdary 
element mesh generation by trans finite mapping method is 
implemented in the form of discrete transfinite mappings. The 
conventional finite element data structure can not provide the 
program the ability to respond to an analyst's request in a 
reasonable amount of time. Application of an existing data 
structure, the winged-edge (W-E) , to BEM-FEM coupled analysis 
system is introduced which can provide the required perf ojrmance , 
The W-E data structure and related concepts are described with 
particular emphasis on its implementation in a combined BEM-FEiyi 
fracture mechanics analysis system. 

A FRActure Propagation code (FRAPCOl , which incorporates 
all the above mentioned concepts, i.e. fracture mechanics, BEM, 
FEM, automatic mesh generation, mesh and attribute editing, 
data base design, computer -graphics etc. is developed and this 
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integration aspect is described. A nurnber of example problems 
are solved to show the capability of FRAPCO and to study 
the validity and desirability of couplina of FEM and BEM for 
application to interactive crack propagation analysis from 
accuracy, time and storage complexities point of view. 



CHAPTER 1 


INTRODUCTION 


1. 1 Fracture Mechanics and Discrete Crack Propagation 

At the most basic level the essential feature of 
fracture phenomenon is the rupture of interatomic bonds 
in the material. In between the phenomenon involves 
nucleation, growth and coalescence of micro and macro 
voids on cracks in the material. At the other end of 
scale, the process is exhibited in the form of separation 
of a full-scale structural part due to propagation of the 
dominant flaw. 

Owing to this highly complex nature of the phenomenon, 
at present there is no consistent single theory dealing 
satisfactorily with all its relevant aspects. The theories 
developed in the past to deal with fracture of solids 
generally treat the subject from one of the three points of 
view, namely, microscopic or atomic, microstructures and 
macroscopic or continuum. 

From engineering point of view, macroscopic theories 
of fracture mechanics are used. These are based on the 
notions of continuum solid mechanics and classical thermo- 
dynamics which provide the quantitative working tools for 
dealing with the fracture of structural materials. It is 
implicitly assumed that material contains some macroscopic 
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flaws which may act as fracture nuclei and that the medium 
is a homogeneous continuum in the sense that the size of 
the microscopic flaws is large in comparison with the 
characteristic micro-structural dimension of the material. 

The problem is then to study the effects of applied load, 
crack geometry, environmental conditions etc. 

Over the years, classical fracture mechanics has grown 
with the objectivity of predicting the initiation of fracture 
and addressing the problem of predicting structural failure 
as the immediate consequence of fracture initiation. But a 
fracture, once initiated, is not necessarily always unstable. 
So, for a propagating crack one must try to find out its tra- 
jectory, wheather the crack stops after a certain distance and 
what must be done to get it going again. In other words, fra- 
cture propagation can be stable in load control sense. Assum- 
ing linear elastic fracture mechanics (LEFMl conditions, this 
implies that stress-intensity is decreasing with increasing 
crack length. Obviously, it is not sufficient just to compute 
stress-intensity factors for an initial crack configuration. 
One has to update stress-intensity factors as crack length 
changes . 

Propagation of multiple cracks is common in realistic 
problems of fracture mechanics. One lias to accommodate 
simultaneous propagation of many cracks. Thus with each 
growth increment of a given crack, a new boundary value 
problem is generated. Displacement and traction boundary 
conditions may change, stress fields are altered and loading 
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may change in direction and intensity. As a consequence, 
propagation of one crack may cause initiation of another 
crack and cracks may influence each other's stability and 
trajectory. 


With possibility of so many factors interacting, 
one must dispense with a simplification often applied to 
fracture mechanics, that of self-similar propagation. 
Curvilinear (mixed-mode) crack propagation is commin in 
real life problems of fracture mechanics. Moreover, if 
one accepts the mixed-mode stress intensity factors, for 
example, and can be present along a crack length, 

then theoretically, mixed mode fracture initiation can occur 
even if ± where is the critical stress inten- 

sity factor, a material constant. Consequently, one must 
incorporate an interactive theory which accurately predicts 
the critical mixture of stress intensity factors which will 
cause the next increment of propagation. 


1.2 Interactive Computer Simulation of Crack Propagation 

Because of complexity and multitudes of difficulties 
in modelling crack propagation, it is currently impossible 
to create a general 'problem in, solution out’ or 'expert 
system' computer code to model fracture propagation. Rather, 
by increasing interactivity of the code, it is possible to 
do numerical analysis tinder continuous, real-time control 
of the analysis through interactive computer graphics. The 
necessary components of such an interactive code, namely, 

r!a 1 . moeTi rfonor’at-i Vkacso rm anfl 
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interactive computer graphics are introduced in this section. 

1.2.1 Numerical Methods in Fracture Mechanics 

In most of the realistic situations, analytical solutions 
can not be carried out, and numerical solutions has to be 
restored to. At present, one of the most generally used method 
having least limitations, is the finite element method (FEM) anc 
another exciting method which has made profound impact on many 
fields is boundary element method (BEM) . 

Finite element method can be interpreted as a piecewise 
application of variational methods in which both domain and 
boundary of the problem are discretized into number of elements. 
The governing differential equations of the problem are then 
approximated over the region by functions which fully or 
partially satisfy the boundary conditions. 

Another possibility is to use approximating fianctions 
that satisfy the governing differential equations in the domain 
exactly but not the boundary conditions. These techniques 
are called boundary solution method or boundary element methods. 
Both BEM and FEM will be explained in more detail in the subse- 
quent chapters but it is interesting to note that they can be 

t 

interpreted as different weighting residual formulations and, 
thus can be coupled in a single unified code to exploit 
advantages of both. 


In both BEM and FEM, the undermined parameters (e.g. 
displacements) represent the values of the solution at 
finite number of preselected points called nodes, on the 
boundary only in case of BEM and on boundary as well as 
interior of the domain in case of FEM. These values of 
primary variable can then be used in a post-processing 
operation to calculate various fracture mechanics parame- 
ters such as stress intensity factors, strain energy density, 
circumferential stresses etc. at desired points on the 
domain and boundary. These parameters are then used to 
construe the actual state of the cracked body and also to 
determine other important parameters such as direction of 
propagation, amomt of propagation etc. to be used in the next 
iterative step of analysis. Thus a totally interactive 
simulation of quasi-static discrete crack propagation can be 
achieved in which user, in conjunction with the program, 
interacts to control the course of analysis. 

1.2.2 Mesh Generation. Data Base Desioh ahd Interactive 
Computer Graphics 

Discrete crack propagation analysis makes much more 
stringent requirements of the FEM and BEM analyses softwares 
than analyses for which these tools are normally designed. 
There are three reasons for this: 

1. Contemporary numerical fracture mechanics is a 
task of modelling inherently non-continuous 
behaviour (cracks and localizations) with inherently 
continuous numerical techniques (BEM and FEM) . 
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2. Because of (1), the geometric and topologic (mesh) 
description of the problem must change with each 
analysis step. Not only must the geometry and 
topology change but the nature of the change is a 
function of the previous analysis results. This 
implies a close coupling between the analyst's 
front and back-end tools. 

3. The geometric and topological complexity of repre- 
senting an arbitrary fracture surface in an arbitrary 
body is something which is beyond the capabilities of 
most commonly used software tools. 

Reason one, above, dictates that the analysis results 
for fracture problems are much more mesh dependent than 
results of more conventional analyses. In fact, to model 
these processes objectively one must remesh the problem at 
each analysis step to reflect the configuration change due to 
crack propagation. Reason two dictates that there must exist 
one piece of software which contains a hybrid of traditional 
pre-and post-processing capabilities. This piece of software 
extracts the results of one analysis and use it to preprocess 
a problem for the next analysis. The third reason dictates 
that to perform numerical fracture mechanics, one should have 
powerful tools for geometric modelling and automatic finite 
and boundary element mesh generation. In order to achieve 
all these capabilities, one must design the program around a 
data structure that will support a fast query time, dynamic 
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changes in problem size, easy identification of features of 
model such as mesh boundaries and material boundaries. Tra- 
ditional finite and boundary element data structures fail to 
provide this kind of stringent requirements and some special 
data structure which has got these capabilities must be 
restored to. 

A necessary component of an interactive simulation 
environment is computer graphics. A computer graphics 
subsystem appropriate for simulation of complex fracture 
events must Involve mesh displays, interactive inputting by 
pointing to some item on the screen, ability to display 
different pre-selected menu items on screen so as to guide 
the user during the flow of analysis. In other words, it 
must provide fast visual feedback to the user to help him 
perform numerical analysis under continuous, real-time 
control . 

All the concepts discussed so far, namely, fracture 
mechanics, finite and boundary element analysis, mesh 
generation, data base design, remeshing etc. combine in a 
coherent program to produce 'interactive simulation' of 
crack propagation process. By 'interactive simulation' 
here is meant numerical analysis under continuous, real time 
control of the analysis through interactive computer graphics. 
By ' control ' is meant that the pace and the direction of the 
simulation are dictated by the analyst rather than by a time- 
sharing operating system and the program algorithms. The 
philosophy behind the generation of fracture propagation 



8 


systeiTt to be described here is to provide a fracture mechanics 
expert with the tools necessary to predict crack growth rate 
and direction at a given instant in the crack evolution. The 
user can, using his own expertise, decide to allow the crack 
to grow further, according to the prediction of the program 
or to over ride the predictions of the program and decide, 
for instance, to test his own theory to calculate the crack 
growth characteristic. The interactivity of the fracture 
propagation systems provide a flexibility by which the user may 
gain understanding into the fracture propagation process through 
numerical experimentation. 

1.3 Review of Previous Works 

1.3.1 Finite and Boundary Element Crack Propagation Systems 

A number of fracture propagation programs have been 
reported in literature. The early version of the programs 
[l],[2j, [3 J relied heavily on the available analytical 
solutions for calculating stress-intensity factors at the 
crack tip and as such, they were not applicable to arbitrary 
structures containing cracks under generalised loadings. 

Truly generalized fracture propagation systems which 
can handle any geometry and loading conditions stem out of 
a number of special crack tip finite elements used for this 
purpose. First of such type of elements either utilised 
a series solution as displacement fvinctions [4J, [5J or 
elementary polynomials of higher order [6 J, [7j as inter- 
polation fxanctions. But these elements have found limited 
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application in predicting stress intensity factors for 
\in-propagated crack geometry and have not been utilised 
for propagation study because of their difficulty in 
satisfying compatibility condition, utilisation of large 
number of nodes, difficulty in remeshing to reflect crack 
growth etc. 

Another type of element which found some use in 
propagation study is so called hybrid elements. Hybrid 
formulation permit choice of two or more assumed fields of 
behaviour in a single element. In fracture mechanics both 
assumed stress-hybrid formulation [ 8] and assumed displace- 
ment hybrid formulation [9 jhave been applied. These special 
elements could easily be combined with those obtained through 
usual displacement finite element formulation, keeping inter 
-element compatibility intact and also, they can be more 
easily made to handle the crack tip singularity in presence 
of plastic deformation. The SIFs are usually treated as 
unknown variables and obtained directly through solution of 
system equations. Atluri et.al. [10 ] developed OR-FLAW, a 
finite element program for such direct evaluation of SIFs 
for user defined flaws in plates, cylinders and pressure 
vessel nozzle corners. OR-FLAW is applicable to a rather 
limited variety of structures because it does not have a 
generalized preprocessing or mesh-updating capability. Also, 
hybrid formulation is very complicated and though they have 
been studied quite extensively, their implementation has 
been quite limited in general purpose finite element LEFM 


nroarams. 
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Simulation of more realistic crack propagation by 
FEM have been achieved by Ingraffea et.al. fll] , [12 ] . 

In these programs, quarter-point element [62 1 are used 
because it gives acceptable accuracy to model crack tip 
stress and displacement fields, does not require any inter- 
polation scheme to extract SIFs and most importantly, 
could be easily generated in meshing and remeshing automa- 
tically. In [ll] , the program FEFAP (Finite Element Fracture 
Analysis Program) has been described in connection with 
modelling discrete cracking in reinforced concrete. Though 
the program has the capability of modelling general crack 
geometry and multiple cracks along with pre-and post processing 
attachments ,it has low-level interaction capabilities. Later, 
a high level interactive version of the same program, FEFAP-G 
has been reported [12 3 • Still, it s pre- and post-processing 
capabilities are limited because it reliers on traditional 
finite element data structure for this purpose. 

Recently, in [l3 3 Wawrzynek and Ingraffea described 
a versatile 2 -dimensional finite element program, called 
FRANC (FRacture ANalysis Code) . It utilises quarter-point 
elements to model crack tip stress-and displacement fields. 

Very powerful pre-and post-processing facilities have been 
made available in FRANC. It utilises many new concepts, data 
structures, algorithms and is the most versatile of all the 
fracture propagation system described so far. 
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Boundary element method has also been used for crach 
surface modelling. Several strategies have been proposed 
for the analysis of crack problems using BEM. These methods 
include representing the crack as a notch, symmetric crack 
modelling, use of special Green's function and flat crack 
modelling. Representing crack as a notch or replacing crack 
plane with symmetry boundary conditions removes the singularity 
in the algebric system of equations which is obtained when 
upper and lower crack surfaces are modelled in the same 
plane [l4 ] . It is, however, limited to symmetric crack pro- 
blems. Representing the crack as a notch increases the 
modelling error due to the notch opening. The special Green's 
fvinction approach [is], [16] has the advantage that crack 

geometry and crack tip singularities are fully embedded. The 
disadvantage is that some two-dimensional and all-three dimension- 
al problems cannot be formulated using Green's function approach. 
Flat-crack modelling [61 3 represents the displacements along 
the crack surface as the relative displacement between the 
two crack surfaces. This scheme has two critical deficiencies 
as a mathematical model for crack geometries, as pointed out by 
Cruse [61 3 . First, if there is no traction on the boundary 
and only crack surface loading, a non-tmique bomdary integral 
equation is generated. Second, two unknown displacement 
variables, i.e., the relative and total displacements, exist 
along the crack. A possible solution for these problems is 
given by Cruse [ 61]. 
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An alternative strategy for modelling the crack and 
subsequent SIF computations is given by Blandford et.al. [ 18] . 

In this, multidomain discretization is used to model two 
crack surfaces in two separate domains which eliminates the 
problems mentioned above. Perucchio et.al. [19 ] extended 
this multidomain boiandary element crack modelling technique 
to three-dimensional cases. 

Later, boundary element method is used in crack propa- 
gation study. Ingraffea et.al. [20 ] has used boundary element 
to model discrete crack propagation in two dimensions. It 
utilises multi-domain BEM and as such, can handle arbitrary 
geometry and crack shape. Gerstle et.al. [ 12] described a 
three dimensional boundary element program BEM-3D which models 
three-dimensional crack geometry by multidomain BEM. It has 
got extensive graphics attachment for doing various pre- and 
post-processing jobs. But none of the boundary element programs 
described here uses any special data structure support and 
as such they are slow in responding to various user's request 
such as mesh editing and remeshing. 

1.3.2 Coupling of BEM’ and FEM 

In order to exploit the advantages of both BEM and FEM, 
the validity and desirability of coupling of these two methods 
have been studied extensively. Such coupling of BEM and FEM 
for static problems was first proposed by Zienkiewicz [ 2l] and 
Shaw [ 22 1 , In these papers, integral equation was used in 
connection with an energy fianctional. The method gives rise to 
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symmetric 'stiffness' matrices of bovindary element region 
thus making it ideal for incorporating into an existing 
finite element software. 

Brebbia and Gergiou [23] presented a direct method of 
combiningthese two methods and they compared symmetrized 
direct stiffness matric with conventional BEM. Kelly et.al. 
[ 24 ] compared the non-symmetric and symmetric stiffness formu- 
lations for potential problems. A critical study of various 
methods of generating boundary element 'stiffness' matrices 
for solid mechanics application can be found in [25 ],[26] #[27] 

Kishimoto et.al. [ 17 ] have applied a coupled BEM-FEM method 
for elastic-plastic fracture mechanics. A few finite elements 
have been used around the crack tip plastic zone and boxmdary 
elements over rest of the domain, as FEM can handle plasticity 
more easily. Stress-intensity factors have been extracted and 
they correspond quite well with previously published results. 
But Kishimoto et.al. [ 17 ] restricted such application to an 
initial crack geometry and crack propagation study by a coupled 
BEM-FEM has not been attempted, 

1.3.3 Mesh Generation 

Current mesh generators employ one or a combination of 
two of the following techniques: 

Cll Node connection approach 
(2) Grid-based approach 

(31 Smoothing procedures 

(41 Blending functions. 
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In node connection approach, nodes are first generated 
inside the dornain by some random node generation techniques 
[ 28 ] , and the elements are then formed by connecting them 
to form triangles [29] or quadrilateral pOj . 

The grid-based approach arises out of the observation 
that grid looks like a mesh and it can be made into a mesh 
provided that the grid cells along the object botindary can 
be turned into elements. Thacker et.al. [31 }are probably 
the first to publish such a mesh generator. Object can also 
be represented as a CSG [32 ]. Yerry and Shephard [33] 
followed a quad-tree encoding. 

Smoothing procedures [34 1 [35 1 too, are cast in the 
form of mesh generators. However, they are normally used in 
conjunction with one of the other techniques to improve the 
shape of the elements by repositioning the nodal location of 
mesh generated. They are basically base on Laplacian mesh 
generator, which is attributed to Wilson fes]. 

The blending function or mapped element approach 
requires an object be subdivided manually into simple regions, 
each of which consists of three or four sides. A mesh can be 
induced in it by mapping a mesh template of tinit traingle or 
square in the parametric space to the region [37 ] , [38 ] . 

This blending technique can be formulated in terms of the elegai 
projector theory [39] to generalize the concept. 

More extensive bibliography on mesh generation tech- 
niques can be found in [40 ], [67]. 
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1.3.4 Data Base Design 

The interactive finite and boundary element programs 
which are to support real time crack propagation analysis 
must have a data structure which should provide — a fast 
query time, dynamic change in problem size, easy identifi- 
cation features of the model such as mesh boundaries, 
transformation of data into a form that is more efficient 
for the computationally intensive equation-solving portion 
of the analysis. An existing data structure, the winged-edge 
(W-E) , originally designed for a very different purpose can 
be used to provide the above mentioned features. W-E data 
structure was introduced by Baumgart [41] to describe the 
surface topology of polyhedra for computer vision. This data 
structure was later used in boundary representation type of 
solid modeller [42] . Different variations of the original 
form of W-E data structure also exist [43], [44] . 

Weiller [44] has established the validity of such a 
data structure for use in boundary representation type solid 
modeller from topological and graph theoretic point of view. 
Later a combinatorial analysis of the W-E and other such data 
structures has been presented by Woo [45 ] . Mantyla [46] has 
shown how such type of data structures could be used in a 
solid modelling environment in conjunction with the so called 
Euler operators for extensive and fast geometric modelling. 

Application of such a data structure in finite element 
is very new and novel. Wawrzynek et.al. [47 ] describes such 
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application for two dimensional finite element application, 
in which a modified winged-edge data structure has been 
established as a robust data structure which can provide 
the stringent performance requirements of an interactive 
fracture propagation system. 

I* 4 Objectives of Present Work 

Development of a comprehensive fracture propagation 
analysis system utilising FEM and BEM is attempted here. 

The specific objectives of this worh can be listed as: 

i) To study of desirability and validity of applying 
a coupled BEM-FEM analysis to interactive fracture 
propagation analysis so as to exploit the advantages 
of both the methods. 

ii) To develop a computer code utilising a coupling 
strategy which permits assembling of contribution 
from boimdary element region with adjoining finite 
elements so that existing finite element software can 
be utilised with minor or no modifications. 

iii) To make both BEM and FEM or tjieir coupled analysis 
formulation available in a single code so that user 
can use any of them at his own discretion thus making 
the debate regarding which method is better for a 
particular problem superfluous. 
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iv) To utilise existing data structure and mesh generation 
techniques so that a unified approach could be adopted 
for interactive generation and modification of both 
boundary and finite element mesh and problem attributes 
in real-time environment. 

v) To study storage and time complexity requirements of 

a coupled BEM-FEM analysis system that is built around 
an edge based data structure and compare it with finite 
element method. 

vi) To integrate the concepts of fracture mechanics, finite 
and boundary element analysis, mesh generation ,data 

base design, remeshing etc. into a single coherent computer 
code to provide interactive fracture propagation analysis 
tools at user's disposal so that he can do numerical 
experimentation with any of his own theory or concept 
for accepting or rejecting it. 


CHAPTER 2 


CRACK PROPAGATION CRITERION 


2 . 1 Introduction 

The determination of fracture initiation from an 
existing flaw requires knowledge of stress intensity 
around the crack tip, determined analytically or numeri- 
cally as fmctions of geometry and load, and the appro- 
priate fracture toughness, a material state property, 
determined experimentally. These parameters combine in a 
theoretical Mixed-mode fracture initiation fxinction analo- 
gous to a multi-axial stress state yield f\mctions of 
plasticity. A number of theories such as energy release 
rate, maximum circumferential stress etc. are proposed, 
though the need of a criterion that can predict failure 
in a simple and unified manner has not been fulfilled 
[48]. 

Recently, Sih [ 49] , [ 50] has proposed a theory of 

fracture based on the field strength of the local strain 
energy density in an element of material ahead of the 
crack or notch tip. The most important feature of this 
criterion is that it is free from physical constrictions 
and sufficiently general for treating engineering problems 
encountered in practice. In this chapter, the basic concepts 
of linear elastic fracture mechanics are described in 
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section 2.2 and then strain energy density criterion for 
incremental crack growth is outlined in section 2.3. 


2 . 2 Linear Elastic Fracture Mechanics 


Linear elastic fracture mechanics analysis neglects 
the inelasticity at the crack tip such as microcracking 
and plasticity, which are considered to be restricted as 
to a sufficiently small volume such that they can be negle- 
cted. 


Considering a crack problem in an infinite domain, 
Williams [si ] has shown that the stresses for the traction 
free crack may be written in terms of an infinite series 
with respect to the polar coordinates r and 0 . Substitution 
of Irwins [52 ] definition of stress intensity factors into 
Williams eigenfunction expressions yields (Fig. 2.1) , 
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in which K^. and are Mode I and Mode II stress intensity 

factors (SIF) respectively . 


Integration of eqn. (2.1) using strain-displaceTnent 
and stress-strain relations yields the following relation- 
ship for displacements relative to crack tip for r << a 
(a = crack length) [53] , 
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in which G = shear modulus, k = (3 - 4v))for plane strain 
andk.= (3 -v)/(l +v ) for plane strain, v = Poisson's 
ratio. 


2. 3 Strain Energy Density Criterion 


Strain energy density criterion makes use of strain 
energy density function dW/dV which in general can be 
determined at any point from the relation [48 ] , 

e , 


w = ^ “’ij 'ij 


(2.3) 
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where and are the stress and strain components 

respectively. In the linear elastic range thus one gets 

dW/dV = a . . e . . /2 . 

ij ij' 

In strain energy density criterion, the first 
attention is focussed on material elements at a finite 
distance r^ from the point of failure initiation e.g. 
the crack tip (Fig. 2.2). The circular region of radius 
r^ represents the core region in which the continuum model 
fails to describe in detail the state of stress and defor- 
mation. The values of strain energy density dW/dV are 
calculated along the boundary of the core region and atten- 
tion is concentrated on the stationary values of dW/dV. 

Thus it is assumed that the direction of the element that 
initiates fracture or yielding corresponds to the minimum 
or maximum value of the strain energy density fvinction, 

(dW/dVl^._ or CdW/dV)^^^. 

min max 

The critical values of dW/dV for yielding and fracture 
are obviously different. Referring to the true stress-true 


strain curve of the material (Fig. 2.3), the critical value 

of dW/dV for yielding, (dW/dV)^ , is equal to the area under 

this diagram upto the point of yielding, while the critical 

value of dW/dy for fracture, (dW/dV)^. is equal to the area 

mm 

of the diagram upto the point of fracture. For brittle 


materials, the values of (dW/dV)*^ and (dW/dV)^, are almost 

max min 

equal, while for ductile materials always 


greater than CdW/dVl‘ 


This fact combined with the previous 


observation that yielding and fracture are associated with 




G. 2-2 Crack tip core 
variation . 







True stress 



FIG. 2-3 A typical true stress - true strain diagram 
in tension for ductile materiol 
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the maximuTTi and Tninimum value of dW/dV respectively, 
suggests that in ductile materials yielding always 
preceeds fracture initiation. 

Based on the above arguments, the basic hypotheses 
of the strain energy density criterion for fracture and/or 
yielding may be stated as [48] , 

Hypothesis (1) : The location of fracture coincides with 

the location of minimum strain energy density rr\±n' 

and yielding with the maximum strain energy density 
(dW/dV) 

Hypothesis C2) ; Failure by stable fracture on yielding 

occurs when (<3W/dV)„,. or (dW/dV)„^„ reach their respective 

TTJin max ^ 

critical values. 


Hypothesis (3) ; The amount of incremental growth r^,r 2 ,-.--r 
...r^ is governed by 


/M) = fi ^ fl 

''dV^c r^^ r2 




constant. 
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where the quantity S is defined by. 
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There is imstable fracture or yielding when the critical 
ligament size r^, a material constant, is reached. The 
quantity S represents the local energy release for a 
segment of crack growth r. It is to be noted that eqn.(2.5) 
is independent of constitutive relations of the material. 

Considering the variation of strain energy density 

function dW/dV versus distance r along the direction of 

expected crack growth (Fig. 2.2), one observes that crack 

initiation starts for the load for which the dW/dV - r curve 

passes from the point defined by the radius r^ of the core 

region and dW/dV = (dW/dV)^. = (dW/dV) . Thus, in Fig. (2.2), 

mm cr 

the crack starts to propagate for q^^ = q _2 for q = q^ 

the crack does not grow and for q = q^, it grows by an amount 
r = r^. Unstable crack growth leading to global instability 
takes place when r = r^, where critical distance r^ is defined 
from relation (2.4). This is equivalent to S = where 
represents the fracture toughness of the material. 

As pointed out earlier, when inelastic deformation 
is negligibly small everywhere in solid, linear elastic 
fracture mechanics can be used to describe the state of 
stress and strain. Thus the quantity dW/dV can be written 
as [49] , 
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Introducing eqn. (2.1) into eqn. (2.6), the following 
quadratic fonn for the strain energy density is obtained. 


If ' r '^11 ■'l + 2^12 ><1 "ll *22 ''ll> 


(2.7) 


where the coefficients a^j (i,j = 1,2) are given by. 


16G a^^^ = (1 + Cos ©) (k - Cos ©) 


(2.8a 


16G a .^2 = Sin © [2 Cos © - (k-1)] 


(2.8b 


16G a 22 = (1^+1) (1 - Cos ©) + (1+Cos ©) (3Cos © - 1) 


(2.8c 


Substituting eqn. (2.7) into eqn. (2.5), one obtains expre- 
ssion for strain energy density factor S, 


^11 ^^12 ®22 


(2.9) 


Thus one can now express hypothesis (1) of strain energy 
density criterion in njathematical forni. 
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3 © " 
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On substitution of eqn. (2.9) into eqn. (2.10), the final 
form is obtained, 
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[ 2 Cos O - (k-1) ] Sin © Kj + 2 [ 2 Cos 2© - (k-1) 

Cos © ]Kj Kjj + (k-l - 6 Cos ©) Sin © = 0 (2.11c 

r 2 Cos 2© - (^^’l) Cos © ]K^ + 2 [(k-l) Sin © - 4 Sin2©]KjKjj 

+ [ (k-l) Cos © - 6 Cos2©] 0 (2.111 

Relations (2.11a) and (2.11b) represents the general 
formulas of the strain energy density criterion and they are 
applied in conjunction with hypothesis (2) and (3) and 
eqn. (2.4), to predict crack growth characteristics. Algorith- 
mically, at any load step, values of stress intensity factors 
are extracted from stress analysis of the body under consi- 
deration and they are introduced into eqn, (2.11) to obtain 
angle ©, which predicts the direction of propagation with 
respect to local crack axis? the crack starts to grow 
only if strain energy density dW/dV, given by eqn. (2.7), 
along the direction of crack extension at a distance r^ 
from the crack tip reaches the critical level (dW/dV)^ and 
grows by an amo\int r, given by eqn. (2.4). 



CHAPTER 3 


A BRIEF DESCRIPTION OF FEM AND BEM AND THEIR COUPLING 


3 . 1 Introduction 

This chapter briefly outlines the theory of finite 
eleinent and boundary element methods in regard to the 
stress analysis problems of fracture mechanics. Since 
finite element method is now well standardized, only a 
brief discussion is given in section 3.2, while bovindary 
element method has been described in section 3.3 in greater 
detail, because of its recent origin. The chapter is 
concluded with a note on the method employed to couple 
BEM and FEM and advantages gained by such coupling. 

3.2 Finite Element Method 

In this section, the finite element foirmulation 
based on weighted residual technique is presented, basically, 
with reference to elastostatics problems. Even though the 
problem could be tackled in many different ways, the weighted 
residual approach seems to be particularly attractive because 
of its generality and also, because boundary element and 
finite element procedures could be treated as essentially 
two different types of formulation from same family of 
weighted residual technique, thus facilitating their coupling. 
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3.2.1 The Weighted Residual Process 

Consider an elastic body defined by domain n and 
boundary r (Fig. 3.1). Its equilibrium is established by 
analysing a system of forces and moments for a generic point. 
Thus, equilibrium equation is obtained, 


a, . , + 

ij /i 


= 0 


(3.1) 


where is the stress tensor and b^ are the body forces. 
The traction and displacement boxindary conditions are given 
as , 


= o 


. . n . 

11 1 


= P. 


on 


r 


2 


(3.2) 

u = u^ on r 

where n^ is the outward normal. 

The general weighted residual statement can be 
written as [54] 
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where u, and p, are the displacenients and tractions corres- 
K K 

★ ★ 

ponding to the weighting field, that is, p^ = n^. 

FEM, the displacement boundary conditions are satisfied 
exactly, thus making the residual on r ^ portion of the 
boundary zero. So the egn. (3.3) can be written as, 

i '"jic.j ^ V >4 ^ f ‘Pk - Pk' ^ ° 

' (3.4) 


Integrating the first term in eqn. (3.4) by parts, the weak 
statement is obtained, 
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dr = 0 


(3.5) 


Limiting now the choice of weighting functions, so that, 




0 on 


\ = -< on Tj 


(3.6) 


and putting egn. (3.6) into egn, (3.5) , the final form of weak 
statement is obtained, 

“jk Mk "k ( Pk “k 

0 n ^2 

(3.7) 
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The independent variable u is approxiTnated by a finite 
series (writing in matrix form) , 


{u } 


M 

I 

m=l 




(3.8) 


where [ are the shape, trial or basis functions, 

is a vectors whose components are the approximations to 

the displacements at node m, and M is the total number of nodes 

The strains at any point are obtained, by applying 
an operator L on {u}. 


{e> = 


M 

1 

m=l 


M 

[N 3{u } = Z [b ] {u } (3.9) 
m** m .'■mm 

m=l 


and consequently the stresses are, 


M 


{a} = [d] {c>= [d32[B 1 {5^1 


(3.10) 


Tn=l 


where [b 3 is the strain- displacement matrix, and [d] is 
the stress-strain matrix. The expressions for l / [b],[d] 
can be found in [54 ] . Now using the shape functions as the 
weighting functions (Galerkin method) , and substituting -for 
stresses and strains, in terms of displacements from eqn. (3,9) 
and eqn, (3,10), eqn. C3.7) can be written as, in matrix form 
C<3ropping the body force term) , 


M 

Z 

m=l 


[ / 


[d][b^1 dn]{u^} = [ Nj] {p} dr 


(3.11) 
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which is a set of algebric equations of the form, 

[K] {u} = {f} (3.12) 

where, the matrices [k] and {f } are obtained by summing 
the contribution from the individual element matrices, whose 
components are defined by 

{f,} = ■'■-e CN®]{p>dr 

^ 2 1 

The presentation so far has been completely general 
and, as such, can be applied to any element form, from a 
wide variety of elements now available. For two dimensional 
stress analysis, 8-noded isoparametric element is very useful 
because it can be used to model geometry with curved boun- 
daries and it gives satisfactory degree of accuracy in the 
solution. This element, details of which can be fovind in [54] 
is adopted in the present application. 

The expression (3.12) represents the overall equilibri 
equations , where matrix [K] is the overall stiffness matrix 
or global stifsness matrix and {f} is the global load vector. 
The matrix [k] is a singular matrix as given in eqn.(3.12) 
and can be solved only after incorporating the specified 
displacement boundary conditions. 


= / r B® ][ D ][ B® ] dfi 

e 

^ (3.13) 
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The solution of eqn. (3.12) constitute a major part 
of the required computer time and has a great bearing on 
the cost of analysis. This has provided a great motivation 
for the development of various techniques. One particular 
technique called frontal solution method [ 55] , which intro- 
duces a natural substructuring or partitioning to make the 
maximum use of uncoupling effect between the nodal degrees of 
freedom, has become very popular. Requirement of very small 
core storage and possibili ty of natural node numbering of 
added degrees of freedom (such as due to propagation of 
crack) without any node renumbering scheme, makes it an 
attractive choice for the present application. Particularly, 
its ability to solve for many load conditions without recom- 
puting and reassembling the global stiffness matrix is very 
useful in crack extension study, where, several load incre- 
ments may be required to drive a crack to the onset of 
propagation. 

The solution of system of eqn. (3.12) gives nodal 
displacements. The element stresses are then calculated 
at the element Gauss-points from eqn. C3.10). 

For analysis of bodies containing crack, the 
standard elements such as an eight-noded isoparametric 
element, fail to give accurate result as they do not contain 
the dominant ^ r term in their displacement inter- 

polation function. Such/ r variation of displacement and 
l//r stress variation arotmd the crack tip can easily 
tio achieved using standard eight-noded isoparametric elements 
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with 'mid-side' nodes displaced from their original 
position to the quarter point locations [ 62] . 


3 . 3 Boimdarv Element Method 

Boundary element method, as the name suggests, 
needs discretization of the boundary of the domain only, 
because it chooses the trial function in such a manner 
that the differential equations are satisfied through out 
the domain a priori. BEM also falls under the weighted 
residual scheme in broader sense, though, as will be 
seen in this section, special techniques are applied in 
its formulation which make it radically different from 
finite element method. 

3.3.1 Boundary Integral Equation 

Here again, the general weighted residual statement 
(eqn. 3.3) is integrated by parts and the constitutive 
relations, i.e. and = | 

are substituted to give/ 

“ i ^jkli "ii Q ’"k "4 
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Intearating by parts again the first term in eg. (3.14) 
and using the reciprocity principle due to symmetry of 

^ijkl ' 

^ ^jk ^jk ^jk '^jk (3.15) 

n “ 

One obtains 

'“jk,!' V“ *' -^2 Pk Pk \ ■3'' 

+-'‘ri \ pj^ <ar +/r^ \ dr (3.16) 

Taking into consideration that the body forces are unknown 
functions, the second integral on the left hand side of 
eqn. (3.16) does not introduce any unknowns in the domain £2 . 
The first integral, however, presents unknown displacements 
in the £2 domain. The objective of boundary element method 
is to eliminate the integral in the domain - the first 
integral on the left hand side - by proposing weighting 
functions which satisfy the equilibrium equation in £2 . 

In other words, one looks for Green's ftmction or fundamental 
solutions satisfying, 


(3.17) 



36 


in which A(s, x) represents the Dirac Delta function, 5 is 

the singular load-point , and, x is the field point . Here 

X and Z belong to a generalized domain Q , which encompasses 

the domain of interest ^ (Fig. 3.2). This effectively means 
* 

that is assumed to correspond to a positive mit load 
applied at point 5 e in each of the three orthogonal 
directions given by the unit vectors e i.e. 

\ ~ A (C, x) e^ (3.18) 

•k 

Furthermore, as the starred displacements and tractions, u^ 

* 

and p^, represent fundamental solutions, so they can be 
written as 





(g, x)e^. 


C3.19) 


4 = ® 1 - 


* * , 

where u^^ (£, x) and ( 5 , x) represent the displacements 

and tractions in the h direction at point x corresponding 


to a unit point load acting in i direction (e^) applied at 
point 5 . 


Substituting eqn. (3.17) and (3.19) into eqn. (3.16) 
and noting the property of Dirac delta function, i.e., 


/ A ( 5 , x) u. (x) e. dn = u. (?) e. t3.2Q) 
£2 1111 


one obtains, in general. 



*2 



and boundary f 



JL* r with the same elastic property 



FIG. 3*3 Singulor point % on the boundory surround 
by part of a circular arc . 
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U. (5) = / u. . (s,x) p.(x) dr (x) -/ p* . ( 5 ,x) u.(x) dr (x) 

+ / u . . ( c , x) b . (x) dQ (x) (3.21) 

J 

Equation (3.21 ) , obtained by weighted residual technique, 
is known as Somialiana ' s identity fss] for displacement 
and can also be obtained using extended version of Betti's 
work theorem. 

Somigliana's identity is a continuous representa- 
tive of displacements at any point ^£^2 for known 
boundary displacements and tractions. The extension of 
it to a boundary point CeT is obtained by increasing 
the domain around the point xmder consideration (Fig. 3.3) 
by a circular region of positive radius and then applying 
eqn. (3.21) and taking limit when £ 0. This gives the 

final expression of boxmdary integral equation [57], 

(C) Uj (Sl+/ Pj^j(C,x) Uj (x) dr (x) = f u^j (^,x) 

P , (x) dr (x) +/ u. . (C,x) b.(x) dn (x) 

J O J 


(3.22) 


Where the coefficient dependent only on the 

geometry of the boundary and is obtained by, 


Cij(C) =«„ 


lim /_ P^j(g,x) dr Cx) C3,23l 


e-»^0 
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Because of the singular nature of the second integral 

on the left hand side of eqn. (3.22) is to be interpreted 
in the Cauchy Principal value sense. The coefficient 

(5) becomes egual to 5^^/2 if the tangent plane at ? is 
continuous; however, for practical applications, it will 
be seen that together with the corresponding principal 

value can be indirectly calculated by applying eqn. (3.22) 
to represent rigid body movement. 

Neglecting body force term, the eqn. (3.22) can be 
written as, 

C^j(5) Uj (C) + /p!j{C,x) Uj(x) dr(x) =f u*j. (C,x)P^. (x) d^ 
r r 

(3.24) 

Equation (3.24) is the starting equation for the boiindary 
element technique. It is valid for both infinite and semi- 
infinite medium and also for bounded bodies as well [S?]. 


3.3.2 Ftindamental Solutions 

Following the definition of fundamental solutions 
Ceqn. 3.17), several different solutions can be applied 
in the formulation. Here, the f\indamental solution due to 
Lord Kelvin [ 563 is presented which serves for a variety 
of problems. This corresponds to tFig* 3.2) to be an 

•k 

infinite elastic medium and consequently, r is taken to 




40 


infinity. 

The expressions for two-dimensional case are 

given by, 


(£,x) = 


-1 


8tt(1-v) G 
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(3.25) 
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where G = shear modulus, v = Poissons' ratio. Above exnre— 
ssions are valid for plane strain case and expressions 
corresponding to plane stress case is obtained by substituting 
V for V where v = v/(l +v ) . 

Also, r = r (S,x) represents distance between the 
load point C and the field point x are its derivatives are 
taken with reference to the coordinates of x, i.e., 

1/2 

r = (r-r^) 

ri = x^Cx) - Xj. (5) (3.27) 

~ (x) ’’ ^i-^^ 
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A.1 though only bounded bodies are considered so far, 
the validity of above fundamental solutions could easily be 
established in case of infinite domain [ 57] . 


3.3.3 Numerical Implementation 

Instead of attempting closed form solutions of 
eqn. (3.24), which is a difficult task and only attainable 
for simple geometries and boundary conditions, the boundary 
element method employes a numerical approach. 

For discretization of eqn. (3.24), the boundary r 
is approximated by using a series of elements , the Cartesian 
coordinates x of points located inside of which can be 
expressed as, 

X?(n) = ^ X, [ d(e, 1) 3N^(n) (3.28) 

where L = number of nodes of the element e, d{e,l ) = global 
e 

number of node L of element e, i = 1,2 : diamension of 
space, n is the homogeneous natural coordinate associated 
with the element e, and l^(n) are the shape functions [57], 
usually polynomial splines. 

The BEM also assumes that the unknown variables 
(u. and p. at the boundary) are approximated in the form 

X 1 
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k=l 
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k=l 




(3.29) 


where N , N are the shape functions associated with the 

node k of the boundary, so that their support coincides 

with the surface of the element to which the point k 

k k 

belongs. Also, u^ and are the unknowns of the problem, 
i.e., the unknown variables at the nodes of the boundary 
elements. For isoparametric elements ,n ,n Nv®^® identical 
and simply denoted here by N. The expressions of N for 
quadratic elements, which is adopted in the present applica- 
tion are given in Fig. 3.4. 

To obtain the values of new M variables of the 
.problem (M = 4 x K) , it is necessary to use a numerical 
algorithm. For BEM, the most frequently used is the point 
collocation which is adopted in present application and 
discussed here. 

Assuming, that discretization of the boundary 
includes P elements, the eqn. (3.24) can be expressed for 
a certain node C , by inserting eqn. (3.28) and (3.29) into 


eqn. (3.24), as 
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C. .(S) u.(C) + I I u [d(e,l)] j P*.[c,x(n) ] 

J e=l 1=1 ■* e 
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N (n) J ( n) dn 

(3.30) 

P Lg 

= 2 I PJd(eA)]/ u* J c,x(ti)]n (T,)J(n) dn 

e=l 1=1 ^ e 
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where J(n) is the Jacobian of the application of eqn.(3.28) for 
each element e. By grouping some of the terms, eqn. (3.30) can 
be expressed as. 
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where. 
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By writing eqn. (3.31) for |= 1,2,...K and collecting 
the different terms, the system of linear equations is obtained. 


[B ] f u> = [a] {d} (3.32) 

v\.hich are 2xK equations with 4xK variables. 

As pointed out in section 3.3.3., the coefficients 
need not be calculated explicitly and can easily be obtained by 
rigid body condition. Assuming tinit rigid-body displacements in 
any one direction for a bounded body, eqn. (3.32) can be written 
as, 


[b 3 { r^}= { 0 } 


(3.33) 


where is a vector defining a unit rigid displacement in 

the direction 1, Hence the diagonal terms of [ B ] , i.e. right 
hand side of eqn. (3.31c) are simply, 
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^ z 5 j f — ^ \ 1 3 el 
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with 


1, C = dCe,l) 


'£ 0, 5 d(e,i) 


C3.34b) 
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For infinite or setniinfinite body, however, a further 
term '•mst be added [ 60l , 


P 

z 

e=l 


L 

e 

I 

1=1 


'^ 5 , d (e • 1 ) 


B 


5 

i jel 


1 


P 

Z 



e=l 1=1 
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(3.34c) 


Once the characteristics of the program and initial 
data have been defined, one can then proceed to compute matrices 
[a],[b] in eqn. (3.32). The introduction of boundarv conditions 
(two at each node) , together with an appropriate assemblage 
algorithm, allow one to obtain a system of linear algebric 
equations , 


[ M ] {x> = {F} (3.35) 

where (x) is a vector in which all the unknowns of the problem 
are included - In BEM the trial functions are not locally 
based like FEM and so the matrix [m] is fully copulated and 
asymmetric. puts severe restriction on the problem size 

and largely dilutes the advantage gained by avoiding the dis- 
cretization of domain and thereby reducing number of variables. 
This has res amany techniques to solve this type of 

problem [ 593 * However, if one ultimately treats the 

BEM region stS one o£ the finite elements by a proper coupling 
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strategy, one can use the standard finite element equation 
solution technique such as a frontal solver. 

In order to get the syetm of eqn. (3.35), it is 
necessary to compute a set of coefficients whose expressions 
are given in eqn. (3.31) . This entails integrating some 
functions which depend on the coordinates of the collocation 
and integration points. This is extremely difficult business 
in BEM and is not as straight forward as in FEM. Generally 
Gaussian integration scheme is used to compute the coeffi- 
cients, though problem arises due to singularity of the 
integrals because of the basic singular nature of the fundamen- 
tal solutions. Although this singularity is saved theoretically 
by considering the integrals in its Cauchy principal value sense, 
it is obvious that a standard Gaussian quadrature process could 
provide inaccurate results when the distance between collocation 
point C and the integration point x(n) on the element e, over 
which the integration is carried out, tends to zero, and, more 
over, when C belong to e, so that the functions to be integrated 
are not bounded. 

Though various sophisticated techniques could be 
used to circumvent this problem such as analytic and semi- 
analytic integration, adaptive numerical integration etc. [ 60 ] 
one simple technique [60] is to use special numerical integration 
scheme for singular element and standard Gaussian quadrature 
for the rest of the elements with increasing number of integratioi 
points as the distance between collocation point and integration 


point reduces. 
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For elements close to the collocation point, a 
rational approach is taken as given below: 

S 1 1.5 0 = 12 

1.5 S 1 5.5 0 = 8 (3.36) 

5.5< S ->0=6 

in which S is the distance between the centre of the element 
and singular node 5 divided by the distance between the 1st and 
3rd node of the element and 0 is the number of Gauss points. This 
criterion was found to be accurate enough without losing compu- 
tational efficiency, eventhough it is not acceptable for severely 
distorted elements. 

For the special case of the singular node being 
coincident with one of the element nodes, the coefficients of 
matrix [B] remains non-singular except for the 2x2 diagonal 
submatrix which corresponds to the singular node. This is 
because, even though there is a singularity of type (l/r) in 
p*^. , the shape functions N(n) are of order r and hence their 
product tend to a finite limit as r 0. For determining 
the singular leading diagonal submatrix, rigid body movement 
Ceqn. 3.34), is used obviating the need to go for complicated 
schemes to compute coefficient (Cl and the principal value 
integral . 
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The calculation of coefficients of niatrix [A]is ^nore 
difficult for singular elements as singularity of type ln(l/r) 
is involved. The terms which do not have this kind of 
singularity are calculated with standard Gaussian quadrature 
scheme with increasing Gauss points according to scheme given 
by eqn. (3.36). The terms involving ln(l/r) singularity are 
transformed by writing, 

ln(l/r) = In (a/r) + ln(l/a), (3.37a) 

so that the integral containing ln(l/r) may be written as, 

11 1 

/ In ( — — )dn = / In [ ? -] diT+ b / In ( — )dot 
jl^ 3r V 3r V Q ot 

(3.37b) 

where b is the Jacobian arising out of the linear transformation 
n to a. In eqn. (3. 37b), a is chosen such that a 0 as 
r(Ti) 0 thus making the 1st integral on the right hand side 
regular. The second integral on right hand side has got integrable 
singularity and can be integrated using logarithmic Gaussian quad- 
rature formula [ 58] . 

3.3.4 Stresses and Displacements Inside the Body 

Once the values of displacements and tractions are 
known on the boundary, one can calculate the displacements at 
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any internal point directly from Somigliana's identity, egn.(3.21). 
For calculating stresses, the eqn. (3.21) is dif ferentiate«f at 
the internal points, to obtain, 

°ijCC) f Pj^(x) dr(x) - / p*j^(C,x) u^^Cx) dr(x) 

^ r 

(3. 38) 


where, for the Kelvin fundamental solution, the new tensors 
are [ 57], 


X) = ^ 


— {{l-2y) (r . 6. . + r , 5., -'r 6. . ) 

iT(l-v)r 


2r^i > 


(3.39a) 


- 2-nil-v)r^ (1-2v) 6.^ r^3^ ^j‘-jk ^i 


(^41, r + r .) 


+ 2v(n. r 4 r - + n. r . r •,.) + (l-2v) (2ni^ r . r . 

IfJ/K J,1,K K,1,J 


+ - C1-4V) } 


(3.39b) 
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3.3.5 BEM for Cracked “Bodies 

If a boundary element method is applied directly to 
a non-symmetrical crack, the resulting boundary equations have 
a singular matrix [6l] , since the method can not distinguish 

between two surfaces in the same plane. Several strategies 
have been proposed to overcome this difficulty such as repre- 
senting the crack as notch, symmetric crack modelling, use of 
special Green's function and flat crack modelling. All of 
them suffer from some inherent drawbacks which prohibit their 
use in general mixed mode cases. An extensive discussion 
regarding is given by Cruse [61 ]. 

An alternative strategy has been proposed by Blanford 
et - al. [ 18] in which multidomain discretization using isopara- 
metric quadratic boundary elements are used (Fig. 3.5) . As the 
two crack surfaces are modelled in two separate domains, so the 
problem of two displacement variables along the crack and non- 
uniqueness of boundary integral equation resulting from crack 
surface loading do not arise. The prescribed tractions and 
unknown displacements along the crack surface are easily repre- 
sented independently in each domain. Further more, the crack 
representation is mathematically correct, i.e. the crack is 
not represented as a notch. 

To properly represent the analytically dominant /r and 
1 ^/jr behaviour exhibited by the displacements and tractions 



N2 a 

N3 a 1/2 ii(n* 1) 

a. Continuous quadratic 
element 



N2 a(l-3/2i2)(U3/2>2) 
N3 = n(9/8n* 3/4) 
b. Fully discontinuous 
quadratic element 


FIG. 3- 4 Different types of quadratic boundary elements 



FIG. 3-5 Multidomain discretization of 
cracked body 
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respectively in the vicinity of crack tip, the classical quadratic 
boundary elements has to be transformed. It has been shown in 
finite element literature [62], that the V r displacement variation 
can be represented using the isoparametric quadratic element by 
placing the mid point node at the quarter point, as shown in 
Fig. 3.6f Using the shape functions defined in Fig. 3.6, the 
variation of the displacements and tractions then become^ 


u. 

1 




+ 



/r + 



i =,1,2 


(3.40) 


In finite element applications, the tractions are 
obtained by differentiating the displacement interpolation 
functions, which results in correct l//r singularity in the 
traction variation. However, in BEM, the displacement and 
traction variations are independently represented. The 
inclusion of l//r singularity in traction interpolation is 
accompanied by multiplying the right hand side of eqn. (3.40) 
by /i/r [ 18] 

= (a];+A^ /r + A? r)A7r = B^/Zr + B? + B? /r 


( 3141) 





which corresponds to the correct form; 1 in the above equation 
is the length of the crack tip eleiDent and the expression A /r 
can be written in non-dimensional coordinate form as shown 
in Fig. 3.6. 

The traction singular auarter point element described 
above is extremely useful to model crack tip boundaries. They 
can be used for mixed mode problem of general nature and can 
be utilized to extract SIF's easily. 

Although several strategies could be adopted to arrive 
at the stress intensity factors, such as J-integral, strain 
energy release rate etc, one particular method which is computa- 
tionally efficient, gives accurate results, requires no inter- 
polation and can be applied in general mixed mode problems is 
the displacement correlation method [is] , [63 ]. Also quarter 
point element fits into this scheme directly. 

Denoting the displacements along the crack axis as u 
(crack sliding displacement, CSD) and the displacements normal 
to the crack axis as v (crack opening displacements, COD) , the 
variation of CSD and COD are, using eqn. (2.2) 
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u = U- + (-3u, + 4u_ 
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(3.42) 
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where the pri^ne indicates that the global coordinate nodal 
displacements have been transformed to the crack tip coordinate 
system defined in Fig. 3.7. 

For a symmetric crack problem (Mode I) . u' = = 0 

and eqn. (3.42) reduces to^ 


V' = (4Vg • - )/r/l + (-4v^ + 2v^ ) rA (3.43) 


Equating the coefficients of r in eqn. (3.43) to that of U 2 
in eqn. (2.2) and taking © = 180°, - 
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(3.44) 


where k = 3 - 4v for plane stress and k = (3-v)/(l+v) for 
plane strain. 

The general mixed mode case is shown in Fig. (3.7b) 
and for © = 180°, interaction of and in eqn. (2.2) is 
decoupled, resulting in. 
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(3.45) 


Algorithmically, to get the SIFs, the displacements 
and coordinates of crack face nodes belonging to the quarter 
point elements need to be flagged, and, retrieved, for each crack 
increment solution. These are then transformed to a simple 
subroutine which codes eqn. (3.44) and eqn. C3.45). 



FIG. 3-6 Traction singular quarter point boundary element 



(a) Mode I crock (b) Mixed -’mode crock 


FIG. 3-7 Element geometries for SIF computation 
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3 . 4 Cp-upling of BEM and FEM 

The most commonly used technique to couple the two 
methods .is to transform the boundary element region into an 
equivalent finite element. This has two inherent problems, 
namely, 

i) The stiffness matrix thus formed is nonsymmetric. 
ii) The equilibrium is not identically satisfied. 

In order to over come these difficulties, several methods have 
been proposed [24] ,[25 ]. The energy approach given by 

K^lly et.al. [24 ] is one of the most attractive and conceptually 
simple. 


3.4.1 Energy Approach 

To generate a symmetric system of equations from 
the direct boundary integral procedure, one has to be minimise 
the energy functional of the form, for the boundary element 
region in Fig, 3.8, 


n 



dr - / {p> dr 

^2 


(3.46) 


If the boundary interpolant for {u} and {p}, which are of 


the form 


{u} = [n ]{u} 

{P> = [m ](p} 


(3.47) 




FiG. 3-8 Domain devided into finite element end boundary 
element regions 
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are substituted into eqn. (3,46), u becomes, 

" = p}^ /F mRn] dr{u }-{u!^ p/ [n]^ {p} 

^ r 2 

(3.48) 

The boundary collocation around the boundary element 
region gives direct boundary integral system of equations (3.32) 
in which rigid body condition is implicitly satisfied by means 
of eqn. (3.33), In order to ensure equilibrium condition, 
eqn. (3.32) is supplemented with an equilibrium equation 

/^{P} dr S /j,[M(r)3 dr{p}= {o> (3.49) 

Linking the direct boundary integral eqn. (3.32) with the 
equilibrium eqn. (3.49), one obtains. 


■[A3' 


[b] [O] 

{ p > 

[03^ 

j 

roF [o]. 

{ X} 


(3.50) 


where {Xl takes the role of a Lagrange multiplier and matrix [ o] 
is given by. 


[O] = / [ M(r) 3 d r (3.51) 

The role of the Lagrange multiplier in eqn. (3.50) is therefore 
to introduce a controlled error into each of equations (3,32) 


vrhich can no longer be satisfied identically if the extra 
condition (3.49) is imposed. 

Inverting eqn. (3.50) and partitioning, one can 
express nodal tractions in terms of nodal displacements, 

{p} = [e3{u} (3,52) 

Putting eqn. (3.52) into eqn. (3.48), the functional tt can be 
expressed as. 


IT 


1 

2 


{u}'^f eT; [Mp [N]dr{u} dr{u}'^ 

r r 


(3.53) 


The minimization of it , which requires 3ir/3u = 0 leads to 


[ k'l { u} +{ f'} = {0} 


(3.54) 


where 


and 


[k‘3 = I [e3'^ f [M3'^p^3<3r+ f[[E3'^ /[m 3'^ [N3dr3'^ 

r r 

(3.55) 

(f'} = - / [ n 3'^ {P} dr (3.56) 

^2 


The matrix eqn, C3.54) is of the same form as the 
finite element stiffness equations (3.12), and so, the 
matrices [k*] and {F»} can be assembled into standard finite 
element systems as contributions from a new element. Compatibility 
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between elements is ensured if the boundary interpolants for 
the boundary integral region are identical to the interpolation 
functions of the adjacent finite elements. 

As pointed out in [64] , special care should be 
taken for corner nodes and at nodes having traction disconti- 
nuity, which cannot be taken into account by a classical 
discretization by ordinary continuous elements. One method 
of solving this is to displace the nodes from the element 
corners, i.e. using discontinuous elements [65] where nodal 
values of tractions are chosen inside the element as shown in 
Fig. 3.4. However, such a discretization increases the 
number of traction nodal values and consequently leads to 
a non-square system of equations and, therefore, a non-solvable one 


To over come this problem, the collocation procedure 

is modified such that collocation points coincide with the 

traction nodal value positions. Also, in order to use the 

integration scheme described in section 3,3.3, displacement 

freedom nodes are also made to coincide with this discontinuous 

nodes. Then after computing matrices [B^ and [a ]t[Bl corres- 

* 

pond to displacements { u } - Fig. 3.4) and in order to keep 
compatibility with a possible finite element zone, one must 
return to the previous interpolation scheme i.e. that of 
continuous nodal displacement scheme. This is simply achieved 

if 

by expressing the inside nodal values {u } in terms of . 
end nodal values {u} , 


{u*> = [s] {u} 




C3.57) 
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The main disadvantage of this coupling strategy is that one 
inversion is required (eqn. 3.52). But the other possibilities 
reported in literature [2?] require solution of set of boun- 
dary equations with multiple right hand side, which does not 
produce any significant amount of saving. An error analysis 
of different coupling strategies has been reported in [25], 
where it is shown that the coupling strategy outlined in this 
section gives better results than the other competiting 
techniques . 

3.4.2 Advantages of Couplihq BEM and TEM 

A combined application of both BEM and FEM becomes 
attractive as each method can be applied in the area where it 
is best suited. By such coupling the following advantages 
of BEM are exploited in a linear elastic fracture mechanics 
problem by using boundary elements around the crack tip; 

1. The inherent advantage of BEM in that all approxi- 
mations are confined to the boundary. Thus it gives 
better results in regions of high stress gradient 
such as a crack tip. 

2. The ease of the use of singular elements to reproduce 
the crack tip singularity. Traction singular quarter 
point elements are very easily and efficiently 
generated from ordinary boundary elements and as 

no domain subdivision is required around the 
crack tip, so the solution becomes less mesh dependent. 
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3. The ease of discretizing a moving crack boundary. 
The error associated with crack propagation 
analysis of FEM due to irregular discretization 
around a crack tip is avoided by using boundary 
elements around a crack tip giving a better 
estimation of the crack path. 

4. The need to discretize the boundary only in 
boundary element region. This reduces number of 
nodes and elements in the problem and so consi- 
derable saving in mesh generation, mesh editing, 
remeshing timings is achieved. 

On the other hand, use of finite elements away from the crack 

tip offers following advantages; 

1. Better efficiency of FEM in analysis of finite body 
problems; for such analysis even quadratic boundary 
elements are sometimes less efficient than the 

FEM [66 ] . 

2. Use of localized discretization, i.e. locally based 
functions which produce narrowly banded final 

set of equations. This contrasts with BEM in 
which fully populated unsymmetric matrices are 
generated which largely dilutes the advantage 
gained in BEM by discretizing boundary only. 
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3. Availability of efficient large number of algebric 

equations solution algorithm. Frontal solver is 
utilised in present application which requires 
very little core storage and also it has resolution 
capability in that it can give solutions for many 
load vectors without recalculating and reassembling 
the stiffness matrices again and again. This is 
particularly useful to determine very quickly 
and efficiently the effect of varying load on 
a propagating crack. 



CHAPTER 4 


INTERACTIVE SIMULATION OF CRACK PROPAGATION 


4 . 1 Introduction 


As apparent frotn discussions of previous chapters, 
crack propagation is a highly complicated and geometrically 
non-linear phenomenon. Because of the complexity of the 
modelling of cracks, it is currently impossible to create 
a general "problem in, solution out" or "expert system" 
computer code to model automatic crack propagation. Instead, 
by increasing the interactivity of the code, its flexibility 
is increased so that user learns much more about cracking 
behavior than he or she could with a 'problem in, solution 
out' type of code. 

A necessary component of an interactive simulation 
environment is computer graphics and related concepts. A 
computer graphics subsystem appropriate for simulation of 
complex fracture events must involve more than simple two- 
dimensional plots, charts and mesh displays. There must be 
capabilities for generating and editing mesh, interactively 
changing problem attributes such as loads, displacement 
boundary conditions and an automatic remeshing capability 
which modifies the mesh to reflect the crack extension. 



To achieve all these, the program must use a robust data 
structure supported by intelligent algorithms which perform 
fast queries and manipulate data base. 

This chapter outlines the mesh generation technique 
in section 4.2 and data base design in section 4.3. Lastly, 
the program integration aspect to achieve all the capabilities 
needed for a fracture propagation simulation, namely, the FEM 
and BEM, fracture mechanics and interactive graphics concepts 
is discussed in section 4.4. 


4 . 2 Mesh Generation 

The blending function approach [38 ] , [ 68 ] is currently 
most popular because of its ability to produce well condi- 
tioned meshes for very complex boundaries- One such blending 
fvinction method, called trans finite mapping method is probably 
the most powerful of all mapping methods. It has been adopted 
in the present application because it can be used for generation 
of both finite and boundary element mesh in a unified manner, 

4.2.1 Transflnite Mapping Method 

Let the problem region R be an open, bounded and 

2 

simply connected planer region in E and S, an open unit 

2 

square (0,1) x (0,1) in R . Let there be a one-to-one 
( univalent) mapping U : S <-*-R such that each point (x,y) e R 
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has a unique correspondent (s,t) e S. In other words, 
goal here is to achieve a curvilinear coordinatization of 
R such that each point (x,y) in R may be uniquely referenced 
by its generalized coordinate (s,t). 

A canonical domain (canonical geometry) means a 
standard domain, such as a rectangle, triangle, circle etc. 
for which a natural (curvilinear) coordinatization is known. 

If the given problem domain R is reminiscent of some 
canonical domain, the objective will be to coordinatize R 
with a generalized system analogous to that of the canonical- 
region. To do so, one widely applicable method can be used (de\ 
oped for’Tiathematical description of free-form shapes in two 
and three-dimensions)^ called " blending functions " methods 
[ 69] . In these, a n-dimensional manifold is expressed as 
weighted (" blended ") combination, of manifolds of dimensionality 
n-1 or less. 

Let ^ be a continuous vector valued primitive function 
which maps S invertibly onto R. F is an unknown fxinction 
and one can identify curves in R as being the images under 
F of constant coordinate lines in S. Also, let P be a simple 
projector such that the vector-valued approximation, P [f] 
is a continuous and tini valent map of S onto R. 

On S, consider a cartesian product partion x 
where ? 0 = < S^< ^ ^ ^l"^* • 
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and let { snd be fxinctions satisfying 

the cardinality property. 


^i = 'ik' = 'J1 


(4.1) 


where is the Kronecker delta. 


Further, define the projectors P and P by, 

S 




M 

S F (S^,t) 0^ (s) 
i=0 


(4.2) 


Pj F] = 


N .V 

j F (S, t^) 'i'jCt) 
j=0 


then the product projector (tensor product interpolant) , 


M N 

P,,P^.[1’3 = I S (S. , t.) 0. (s) y . (t) (4.3) 

St i=0 j=0 111 3 

^ M,N 

interpolates F at the (M+1) (N+1) points ((S. ,t.)} 

^ J i=0,j=0 

In contrast, the transfinite interpolant, defined by 
Boolean Sum, 


Pg © ^t f - ^s f Pt t ^ 3 ■ ^s^t [ ^ 3 

(4.4) 

interpolates to ? at all points along the two sets of lines 


{S - Sj[}^_Q and (t - 


N 


The f motions 0j_ and ’Pj are 
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called blending functions . 

If the blending functions are chosen as the Lagrange 
polynomials (of minimal degree) satisfying eqn. (4.1), then 
p, © [ F ] is the transfinite Lagrange interpolant . In 

particular, if M = N = 1, then P © [p] is bi linearly 

blended transfinite Lagrange interpolant, 

Ps ©P^[f]= (1-s) F (0,t) + s F(l,t) + (1-t) F (s,0) 

+ t F (s,l) - (1-s) Cl-t) F (0,0) 

- (1-s) t FC0,1) - s(l-t) F(1,0) -St F(l,l) 

(4.5) 

The images of lines s = constant and t = constant in 
S are curves in R. In particular if boundary of R, 3R is 
identified as being comprised of the four curve segments 
F(0,t), F(l,t), F(s,0), and F(s,l), then 3S 3R for any 
choice of fmctions {0.} and {'?.> which satisfy eqn. (4.1). 

J 

Higher order interpolants may be used to force the co-ordinate 
curves to pass through specified curves on the interion of 
R [68] . 

The general transfinite mapping describes an approximate 

surface or volume which matches a desired or true surface or 

denumerable set of points. It is this property 


volume at a non-' 
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which gives rise to the ter^n transfitiite mapping. This 
property contracts with the isoparametric mappings [ 74 ] 
which match the true surface at only a finite number of 
points. So transfinite mapping can be made to exactly model 
all region boundaries and no geometric error is introduced 
by the mapping. 

Two classes of representation of the boundary of a 
region is possible, namely, continuous form and discrete form. 
The continuous forms represent the position vector to a 
boundary curve as a function of some parametric co-ordinate. 
The discrete representation [ 38] consists of finite lists of 
points located on the curve, with a unique co-ordinate asso- 
ciated with each point in the list. The discrete form of 
representation is completely general and the present precessor 
utilises this concept. 


4.2.2 Implementation 

Discrete transfinite mapping [38] may be used to map 
any mesh in the appropriate primitive polygon to the actual 
region. Thus any mesh topology may be used. However, it is 
convenient to sacrifice topological generality by choosing 
the intersections of the constant coordinate curves for 
use as nodal points in order to minimize input. In bilinear 
mapping, s direction curves are defined by constant t 
coordinate values and vice-versa. Thus a list of n constant 
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t values t^., j = l,n defines n curves running in s direction 
and a list of m constant s values s^ , i=l,TT, defines a set 
of m curves running in t direction. Examination of eqn.(4.5) 
reveals that to evaluate the location of nm intersection 
points of these curves, the boundary curves need to be 
evaluated at 2 (n+m) points: 


{?( 0 , 


tj.), F(l,t^)} i=l,n and {F Csj,0), E(sj,i)} 




For efficient evaluation of all nm intersection points, 
2 (n+m) boundary curve points are evaluated first and stored to 
avoid redundant calculations. 

An apparent drawback of the discrete form of curve 
representation is the large amomt of data required to describe 
a boundary curve in discrete form. This drawback is eliminated 
by providing a package of curve generator that automatically 
create discrete curve descriptions based on various continuous 
models. Each curve is assigned a unique curve number and 
list of discrete points representing it is stored in separate 
arrays . This allows user input to requirement to be reduced 
to a minimum. 

When highly distorted boundary curves are used to 
define aregion, a problem of overspill may result [38 ] . 

This problem occurs when finite element mesh lines cross 
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outside the region boundaries. Higher order ^napping may be 
used to force the mesh lines to pass through some predefined 
constraint curves on the interior of the region [ 69] . In 

the present preprocessor* this problem is overcome by breaking 
the highly distorted region down into two or more regions with 
more regular geometry. Appropriate book keeings are maintained 
to avoid generation of duplicate nodes along the curves shared 
by more than one region. Regions are first reordered so that 
regions sharing a common point may be numbered in a continuous 
sequence. If there is still any duplicate node, the node 
numbers are compressed and the node numbers associated with 
individual elements are corrected. 

In order to apply a uniform algorithm for generating 
both boundary and finite element mesh, regions which are to 
be discretized by boundary elements are treated as a 'super' 
finite element, i.e. a big finite element with more than 
eight nodes. These regions are flagged and while transfinite 
mapping is applied to generate nodes, the internal (domain) 
node generation are suppressed for such regions, so that 
only botmdary is discretized. A finite element type node- 
element connectivity for the super finite element is then 
generated and this is later translated into a boxmdary element 
type node-element connectivity by appropriate coding. 
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4 , 3 Data Base Design 

The program data base is the Kernel which ties the 
rest of the program together and it has a profound impact 
on the design of the rest of the program. The data base 
used in the present implementation is designed around a 
winged edge data structure, introduced by Baumgart [41] to 
describe the surface topology of polyhedra. 


4.3.1 Definitions 

At this point a few definitions are necessary. 

Topological information is the adjacency description 
of a polyhedron. In finite element sense, this would be 
information such as which nodes are adjacent to a given 
element or which elements are adjacent to a given node. 

geometric Information , in contrast, is information 
about actual point locations (l.e. the nodal coordinates 
in a finite element mesh) . 

A vertex is a unique point in space; in planer case, a 
point in the plane in which the mesh construction will take 
place. In the discussion to follow, a vertex is isomorphic to a 
node in finite element sense and the terms will be used 
interchangebly . 
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An edge or arc is a set of two vertices. The edges 
discussed here will be directed edges, i.e. edge points 
from one vertex to the other. The finite element analogy 
to the edge is a portion of an element boundary that connects 
two adjacent elements. 

A graph is a set of vertices and distinct edges 
which utilize the vertices. One can view finite element 
mesh as a graph. A graph that allow self loops (i.e. graph 
configuration in which an edge joins a vertex to itself) 
and multiple edges (or multigraph i.e. a graph configuration 
where multiple edges are allowed to join the same two vertices) 
are called pseudographs . When graph is referred to here, 
it is not referring to pseudo graphs. 

A graph can be embedded (or mapped) onto a surface 
if it is drawn on the surface so that no two edges intersect, 
except at their incident vertices. 

Faces are the regions of a surface defined by a 
graph embedded on a surface. Each face is a connected 
component of the set obtained by subtracting the vertices and 
edges of the embedded graph from the surface. The boundary 
of a face consists of those edges and vertices of the 
embedded graph whose every part touches upon the face. A 
simply connected face has a single, connected boundary. 

A multiply connected face has boundary that consists of 


two or more disconnected components, as in a face with a 
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hole in it. Here, only sirnply connected faces are referred 
to. 

Finite element mesh is a collection of faces, i.e. 
all the elements in a finite element mesh will be faces. 
There is also one special face that arises from the process 
of mapping a polyhedron on to a sphere (to be described 
in the next section) . The face that covers all area exterior 
to the mesh must be stored explicitly with the rest of the 
faces. For such a collection to define a valid mesh, 
the following must be satisfied [70 3 : 

(1) Faces of the collection may intersect only at common 
edges or vertices (if at all) 

(2) Each of the edges is shared by exactly two faces 

(3) Faces around each vertex form a single circuit, i.e, 
they can be arranged in cycle, such that each pair of 
consecutive faces in the cycle meet at an edge adjacent 
to the vertex. 

4.3.2 Winced Edge Data Structure Theory 

At a first glance, it is not obvious that a data 
structure developed to store the topology of a polyhedra 
would be useful .for representing a two-dimensional finite 
element mesh. However, the topology of a polyhedra can 
always be mapped onto the surface of a sphere (Fig. 4.1). 



The mapping can be performed in such a way that one of the 
polygons covers a very large percentage of sphere's surface. 
Attention can be restricted to the portion of the surface 
that contains all other polygons. This essentially is the 
topology of the polyhedron mapped on to a plane, provided 
the curvature of the sphere's surface is sufficiently small. 
Finite element analysis is concerned with this planer topology; 
the polyhedron that is represented by it is not important [47 ]. 

The winged edge polyhedron representation contains 
three sets of records: one set each for vertices, edges and 
faces. Each record contains information associated with the 
topological entity it represents and also contains pointers 
to other records in the data base (Fig. 4.3) . It is by means 
of these pointers to other records that the adjacency infor- 
mation inherent in the data structure can be extracted. The 
winged edge data structure and similar representation are 
called edge based data structures [44] , because the edge 
records are used as keys to access other type of records. 

The edge records are used for this function because the 
quantity of adjacency information associated with an edge is 
invariant from edge to edge. That is, each edge is always 
adjacent to exactly two faces and exactly two vertices, whereas 
a face or a vertex can be adjacent to any number of edges 
(Fig. 4.2) . 

As seen in Fig. 4.2 , the winged edge structure 

maintains the adjacency information with pointers to two faces, 




FIG. 4-1 Mapping of polyhedron (in this case a cube) 
on to the surface of a sphere 
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FIG. 4- 2 The winged-edge data structure 
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"the two vertices and some of the edges adjacent to the reference 
edge. This last set of adjacencies is divided into two . sections 
each associated with the use of one of the sides of the reference 
edges in the circuit of edges around a face. The only edge 
adjacencies represented, therefore, are the four edges that 
directly follow or preceds the reference edge in edge cycles 
surrounding the two faces adjacent to the reference edge. Thus, 
the ewe and cewe field names used here refer to their use in deter 
mining the cycle of edges surrounding a face, as viewed from 
outside the plane looking towards the surface. Finally, Neptr 
and Peptr are the forward and backward pointers for a doubly 
linked list that contains all the edge records (Fia.4.3a). 

Figure (4.3b) shows the structure for a vertex record. 

The vertex records not only contain adjacency information but 
they also contain geometric information (i.e. nodal coordinates) 
and information used in the finite element processing. The 
Eptr field points to one of the edges to which the vertex is 
adjacent. The Neptr and Peptr fields are forward and backward 
pointers to a doubly linked vertex list similar to the edge 
list. The vertex record also contains global finite element 
node number associated with the vertex. This information is 
used while assembling the finite element equations to interpret 
the structural response information. 

The structure for a face record is shown in Fig. (4.3c) . 
The record contains the pointers to the face list and a pointer 
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7Q 


to one of its adjacent edges. In addition to the adjacency 
information, the record also contains bounding area information 
for the face. This information is used to speed up the process 
of determining if a given point is or is not within a face. 

If the point is not within the bounding area of the face, the 
face can be rejected and no further processing is necessary. 

The face record also contains a face-type field. This field 
tells weather the face is being used to represent a finite 
element, and if so, what type. The two remaining fields are 
only defined if the face is representing a finite element. The 
face-number field contains an index into data structure that 
maintains all the stiffness matrices of all the elements in 
the mesh. The material field identifies the elements' material 
type. 


4,3.3 Euler Operators 

The main problem of boundary representation type solid 
modeller, for which W-E data structure was first introduced, 
is their complexity; models of practically interesting objects 
consists of tens or hundreds of faces, edges and vertices. 
Consequently, their construction is an elaborate and error- 
prone task. 

Euler operators [423 were first introduced in order 
to simplify the description of boundary representation. Their 



key idea is that the construction of each meaningful 
boundary representation can take place in a step wise 
fashion using a small set of data structure constructors 
which effectively hide to details of representation 
used. 

Consider a simple polyhedron with V vertices, E 
edges and F faces. Then, according to Euler formula, 

V - E +F = 2 C4.6) 

Euler operators derive their name from this formula. 
Following Braid et.al. f 42], eqn. (4.6) can be interpreted 
as the equation of a two-dimensional hyperplane in the 
three dimensional lattice axes < V,E,F> . Obviously, the 
two-dimensional plane can be spanned by a collection of 
two-linearly independent vectors of the hyperplane satisfying 
eqn. (4.6). Euler operators then may be viewed as such 
vectors, which add or remove vertices, edges, faces as deter- 
mined by the components of the corresponding vector of the 
plane given by eqn. (4.6). List of operators chosen for 
the present implementation are given in Table 4.1. 

The semanties of Euler operators are best represented 
by plane models. Only informal descriptions of each operators 
are given with reference to Fig. 4,4? rigorous formal 


Table 4.1 ; Euler Operators 
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definitions can be found in [70 ]. 

Operator : MSFV forms the initialization operator 

set of vertices, edges and faces. It creates a new 
graph instance with just one node and face. MS^v is always 
applicable. 

The corresponding inverse operator KS'F'v is applicable 
to a graph component consisting of just one face and one node. 
It removes the face and the node. 

MEV Operator ! The MEV operator subdivides the cycle 
of edges of a node, thereby creating a new edge and a new node. 
It basically acts as a "draw-an-edge" operator. 

The inverse operator KEV is applicable to an edge 
occurring twice in a region and having no attached edges at 
one of its nodes. It removes one edge and one node. 

MEF Operator ; The MEF operator divides a face into 
two faces by adding a new edge between two nodes. The new 
edge will occur once in both new faces in reverse orienta- 
tions . 

The inverse operator KEF is applicable to an edge 
occurring in two distinct regions. By KEF, the two faces 
are connected by removing the edge to form a new face 
bounded by the remaining edges and nodes of the original 
regions. Note the duality between MEF and MEV, 



MSFV 

Null 

KSFV 

(a) 





Fig.4*A Euler operators 



SEMV Operator ! The SEMV operator splits an edge to 
introduce a new node and a new edge thereby adding a new 
node and a new edge to the face cycles. The new edge and 
node will occur in both the faces to which the original edge 
belongs . 

The inverse operator JEKV removes an edge and a 
node by joining two nodes by a new edge. It changes the node 
and edge cycles of two faces to which the edges belong. 

4.3.4 Adjacency Relationships 

Adjacency topology concerns the physical adjacencies 
of three primitive topological elements (vertices, faces 
and edges) on the graph. To facilitate discussion regarding 
this following symbols are used: 


V. a vertex 

1 

an edge 

F. a face 

1 

{V^} a set of vertices 

(e.) a set of edges 

1 

{F^} a set of faces 

V total number of vertices 

E total number of edges 

F total number of faces 


Since a three element topology ie aieonseed here 

distinct combination of element types rn the 
there are nine distincu 
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adjacency relationships; each is called an adjacency rela- 
tionsh'.r class or topological query (TO)- These nine TOs 
are in Table 4.2. Note the duality between a vertex 

and a face. In particular, note that EE^ wings of an edge 
are self dual. 

h ooundary data structure can be thought of as a set 
of ad-'acencv relationships among the three basic topological 
entitities [73] . Let a relation be denoted by 


(4.7) 


where X,y can be a set of vertices {V . } , edges , 

faces . A relation {E.}->{V.> or simply E+V, for 

example, stores two vertices for each of the eSjes. Hence, 
given an edge E., its associated vertices can be 

retrieved or updated. All such hlne adjacency relationships 
or TOS can be represented as a graph with three nodes and 
nine arcs, as shown In Mg- 4.5a. Figure 4.5b shows the 
schena for wlnged-edge data structure, with an edge pointing 

two faces, and four of the possible 
to its two vertices, two 

^ tr-K-n^ a face or a vertex points to one of their 

many edges ^ while a 

«= Note that three TOs e.g. E V, E F and 
many edges • Note 

a. -ori in W-E data structure. Some TQs 
E E are fully stored xn w a aa 

tionally e.g. instead of storing all the 


are stored fraction 
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(a) Nine relations and three 
entities 
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(b) Winged edge data structure 

FIG. 4- 5 Schematic for edge based data 
structure 



EF^ edaes for a face F^, a face pointing to only one edge 
is stored. Similarly, V -»■ E is also stored fractionally. 
These two and other four TQs which are only partially 
stored or are not stored at all, are to be accessed by 
intelligent algorithms which make use of information access- 
ibility available in the W-E data structure. Such algorithms 
are discussed in detail in [44 ]. 


4.3.5 Storage and Time Complexity of W-E Data Structure 

The storage complexity of a data structure is measured 
using counting formulas, using which one can assign each 
relation a storage cost in terms of the total number of 
edges E in the graph. Therefore, a set of adjacency relations 
with a storage cost represents a static view of the data 
structure. Also by the basic queries for accessing, one can 
describe a relation that is not stored directly in a given 
data structure and that can be expressed as a procedure in 
terms of relations that are stored. The time complexit y; 
of a data structure is measured using this set of primitive 
queries. Hence, here one gets a dynamic view of a given 
data structure by the way it is accessed indirectly or 

reversely. 

The storage complexity of a relation X Y can be 

[■45] , 


measured by taking the sum 



X 

z 

i=l 


(4.8) 
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YX. 
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Where X and Y can be V, E or F, and i is suT^med over all X. 
Implemented as arrays, the storage for two relations E V 
and E ^ F cost 2E cells each. This is because each edge E^ 
has two vertices vert [ 1 3 and vert [ 2 J as well as two faces 
face [l 3 and face [ 23 . As there are E such edges, the 
total storage is 2E. Hence, 

E E 

I VE. = Z FE. = 2E (4.9) 

i i ^ 

As each edge has four wings, the storage cost for relation 
E -»■ E is clearly 4E, i.e., 

E 

^ EE. = 4E (4.10) 

i ^ 

In W-E data structure, use of fractional relation is 
made, i.e., instead of storing all the EF^ edges for a 
face F^ , a face pointing to only one edge is stored. 

This fractional relation F E gives rise to a storage cost 
of F as it points to only one edge, and there are F of 
them. Similarly, V E costs V in storage. The total 
storage is therefore 8E + (V+F) = 8E + (E + 2) = 9E. 

[ making use Euler formula, eqn. (4.6) 3 • Implication of 
this is less the number of edges in the graph, less is the 



storage requirement. 


The time complexity for the W-E data structure can 
be analyzed as follows. Since the three direct relations 
are E -»■ V, E F and E -*■ E, the three corresponding queries 
T04/ T05 and T06 can be answered in constant time K as 
the arrays allow direct access. For other six queries , local 
'clocking' around a face or around a vertex is required, 
costing EF^^ or EV^ in time respectively. Since E is 
of the same order as EF^ ,one can chose EV^ as the unit [ 723. 
Hence, the total time complexity for the winged -edge is 
0(6EV^ + 3K) . Implication of this is that, time complexity 
is dependent upon a local quantity EV^ instead of a global 
E. In other words, processing time is independent of the 
size of the problem (which is equivalent to total number of 
edges E) , i.e. response time will not suffer as the problem 
size increases. Also, if the local quantity EV^ decreases 
(i.e. less cluster of edges around vertices, in general), 
access time falls linearly with EV^ making data processing 

time less . 


4,3,6 Implementation 

The program here is written in FORTRAN/77 programming 
language. A serious drawback of using FORTRAN is that the 

language does not support record-structured data types, 

«i-nre the winged edge data structure 
the natural way to stor 



information, with FOR^N, the subscripting mechanism is 
chosen for providing the necessary addresses. A set of 
FORTRAN vectors is used to represent the various fields com- 
prising the records in a linked list structure. A single 
record is represented as a cross-section of these vectors 
across the same subscript, and this common subscript serves 
as a pointer to that record. Consequently, the vector names 
serve as the field selectors, while the pointers are used 
as subscripts for various vectors. 

The software support for the data structure is developed 
in a number of distinct layers to ensure modularity (Fig. 4.6). 
The lowest level contains the memory allocation and deallocation 
routines. Control over storage utilization is maintained by 
keeping all the \inused records chained together into available 
lists. Whenever a new record is to be created, memory allocation 
routine remove the first available record from the availability 
list. When a record is deleted, it is added to the avail- 
ability list so that it is reused before any new memory is 
employed by the data base. 

The next level of software contains the routines to create 
and delete the elementary topological entities. These routines 
create and destroy vertices, edges and faces. For these 
processes, the record creation routines are called so that 
memory is allocated for the new record or memory passed onto 
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Euler Operators 


Create/Delete Edge, Vertex, Face 


Memory Allocation/Deallocation 


Data Storage Memory 


Fig. 4.6 : Layered Organisation of Data Base 

Modification Routines. 
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these record creation routines. 

The software level above the routines that deal with 
individual topological elements are those that implement 
Euler operators described in section 4.3.4. These routines 
call the routines which create or delete edge, vertex, face 
to achieve various functions such as making or killing edges, 
vertices, faces etc., keeping topological integrity intact. 

Also they operate on few other fields of certain records e.g. 
material code field of face record. 

In addition to creation and deletion routines, the 
data structure software also supports a number of query 
routines. All the query routines of Table 4.2 are supported. 

In addition, the following additional routines are also 
supported: get the list of all vertices, get the list of all 
faces, GET-FACE. The GET-FACE returns the face that contains 
a given X and Y coordinate pair. This allows the program to 
determine which element a user has selected with the cross- 
hair cursor. 

A uniform approach is taken for applying such data 
structure to both BEM and FEM. Boundary element regions are 
treatea^super’ finite element type faces and they are 
flagged with a special face type which differentiates them 
from the possible adjoining finite elements. All queries 
and Euler operators can then' be applied on such a face without 
any restriction or modification. 



4-4 PrograTn Integration and A Sample Problem 


The key aspect of the present FRActure Propagation 
code (FRAPCO) developed here is the integration of all 
topics discussed so far, namely, fracture mechanics, BEM 
and FEM, mesh generation, data base design, into one coherent 
program. The current version of the program operates on 
supermini computer ND-560/CXA system which consists of 
ND-500 , 32 BIT main processor and ND-100, 16 BIT front 
end processor. The system has virtual memory capabilities 
which obviates the need of using secondary storage units, 
i.e. all the program arrays are held in 'core* and the 
swapping of part of any array to fa$t disk storage is left 
to the efficient system routines. The user interaction and 
display relies on hardware and software facilities of 
Tektronix 4109 (19", 640x480 pixels, 16 colours/4096 hues ) 
and Tektronix 4107 (14", 640x480 pixels, 16 colours/6 4 hues) 
raster graphics terminals. 

An important feature of FRAPCO is the modular nature 
of its basic program design. A common data base is maintained 
containing a complete current description of the structure. 
This data base is accessible through data access routines, 
which is used by all higher level routines for query into the 
data base or to modify the data base. 



The systern is menu driven. A schematic of inter- 
action of the menu pagss is indicated in Fig. 4,7. At 
any point of execution of the program, the user is confronted 
only with the options on the currently displayed menu page. 

This greatly helps the user who is guided by the menu page 
structure which is designed to help user tackle particular 
aspect of finite and boundary element crack propagation 
s^^lysis • The finite element mesh and problem attributes 
change during the flow of the program either automatically 
or manually, under the control of the user, through graphical 
and/or alphanumeric key board inputs. Geometrical information 
can be input by using the cross-hair cursor. Previously defined 
elements of a data base can.be identified by pointing with 
the cross-hair cursor. This frees the analyst from the burden 
of referring to nodes, elements, lines etc. by number through 
keyboard. 

The operation of the system can best be appreciated 
by considering fracture propagation analysis of a practical 
problem. For this purpose, a rectangular plate with a crack 
emanating from a circular hole (Fig. 4.8) has been selected. 

The step by step solution of the example problem illustrates 
how W-E data is used to support some of the functions performed 
by the program FRAPCO also how the program interacts with 
the user for a crack propagation analysis. 




F16. 4-8 Crack emanating from a hole In a plate subject 
to uniaxial tension 


4.4.1 Geonietrv Handling 


'T'he first step in the analysis is the Qeometric des- 
cription of the problem. Geometry handling refers to the 
description and modification of the geometric and topological 
aspect of the problem. Under this come the mesh creation 
and mesh modification routines. The first step in this is 
to specify the data limits i.e. the extent of world window 
which is to be transformed into device (i.e. display terminal) 
coordinates system for displaying the mesh. 

The boundary curves defining the subregion outlines 
are input by the curve generating package invoked by the 
outline menu page, shown in Fig. 4.9. In the present program 
the user can select from the following curve types - straight 
line , circular arc or discrete. The discrete curve option 
is used when mathematical form of the curve is not explicitly 
known. Depending upon the key information (e.g. for straight 
line - starting and end point coordinates , number of sub- 
divisions on the line etc.) provided by the user, the 
program generates and stores discrete form of the outlines. 

If the user .is not satisfied with any particular curve, he 
can delete it through curve delete option by picking the 
curve to be deleted. Fig. 4.9 depicts the complete boundary 
curve description of the present problem. 

The actual element mesh is generated using the display 

page shown In Pig. 4.10. The bowndery curves are aisplayed 
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for selection by the user, who is promised to pick the 
four curves defining each region. The user then chooses 
an element topology (3-noded triangle, 4-noded quadri- 
lateral, 6-noded triangle, 8- noded quadrilateral, 
boundary elements) . However, even though all such types 

of elements could be generated, at present, for numerical 

and 

analysis, only 8-noded quadrilateral finite elemen biquadratic 
boundary elements are implemented. The program automatically 
computes and displays the mesh for approval of the analyst. 
The domain as well as boundary nodes are created for a 
region which is discretized by finite elements and boundary 
nodes only for a region that is discretized by boundary ele- 
ments. A uniform code is used to generate mesh for both 
finite and botindary element regions and no extra effort is 
required from the user to generate a coupled boundary and 
finite element mesh. 

( 

Display of finite element mesh is very efficient task 
with the winged-edge data base being maintained. It is 
desirable that each edge in the mesh be drawn only once. 

If the only information available to the program is the 
element connectivity table, then a check must be made for 
each edge to determine if it has been drawn yet. This 
algorithm will take O(n^) time where n is the number of 
elements in the mesh. With W-E data structure, displaying 
the mesh is a simple ^ since each edge is stored uniquely, 
the program can traverse the linked list of edges embedded 
in the data structure. This process takes 0(m) time where 
m is the number of edges in the mesh. 
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Fig. 4. 11 : Mesh Editing of the Example 
Problem 
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After the .nesh is created, the user may wish to 
perform local mesh modification. This is accomplished by 
using the mesh editor display page shown in Fig. 4.11. 

All the six Euler's operators described in section 4.3.3 
are provided for this purpose. The program performs the 
modification and query tasks through these routines to 
directly work on the winged edge data base. In the 
present example, more boundary element sub-domains are created 

around crack tip in a way to refine the mesh around the 
crack tip. The edges are created by HEV and mEF functions 
to arrive at the final mesh shown in Fig. 4.11. 

4.4.2 Attribute Editing 

After the geometric description of the problem, before 
an analysis can be performed, additional input data such 
as material properties, boundary conditions, loads etc. 
are required. This data may be viewed as a series of 
attributes associated with the nodes and elements contained 
in the geometric description. This inputting is achieved 
through 'Attribute' menu page shown in Fig. 4.12 which 
has options for specifying both point load and distributed 
edge (pressure) load, displacement boundary conditions and 
material properties. 

The process of specifying displacement boundary 
conditions is considerably simplified by assuming several 
nodes on the boundary to have same type and amount of 
displacement boundary conditions. Thus the user is prompted 
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to input the specified displacements and fixity conditions 
to be associated with a particular portion of the boundary. 

The user then indicates the portion of the boundary over 
which these conditions are to apply by picking the starting 
node and end node of the boundary portion. In the present 
example, the bottom edge of the plate is to have all fixed 
displacements in both X- and Y-direction. This 'both fixed' 
code is entered first through keyboard. The portion of the 
boundary over which this condition is to apply is then 
indicated by picking the right-most node (Fig. 4.12) and 
left- most node (Fig. 4.13) of the bottom edge of the plate. 

The program automatically picks all the bo\andary nodes along 
that portion of the boundary and makes all such nodes 'both 
fixed' type. The load boundary conditions are input similarly 
in turn by specifying point load and/or edge load to be applied 
on a group of nodes along the boundary. 


If no boundary information is explicitly available to 

the program, this type of boundary portion identification 

requires an 0(n2) (n = number of elements In the mesh) search 

to identify boundaries. With W-E date-structure, faces in the 

topology that do not represent elements are flagged with a 

special face type. The set of edges adjacent to these faces 

comprise the structural boundary. Once a starting edge or 

node on the boundary Is Cgraphically) identified by the 

analyst, the wing fields in the edge records are used to 

ordered list of edges and nodes along the boundary 


traverse an 




Fig. 4. 12 : Attribute Editing of 
the Examnle Problem 


Fig. 4. 13 ; Attribute Editing of 
the Example Problem 



Fig. 4. 14 : Crack-tin Definition of the 
Example Problem 


104 


spanning fro^n the indicated starting vertex or edge to 
the indicated stopping vertex or edge. 

4.4.3 Crack Tip Definition 

The program FRAPCO uses quarter-point boundary 
elements to model the singular stress field aroiand the 
crack tips. This requires modification of coordinates 
of all nodes adjacent to the crack tip nodes. This is 
achieved by the FRACTURE option in the driver menu page, 
as shown in Fig. 4.14. 

Typically, the user picks the crack tip node by 
pointing to it with the mouse. The program then identifies 
all nodes adjacent to the crack tip node, i.e. the 
'quarter-point' nodes. All the nodes adjacent to these 
quarter point nodes are also identified so that quarter 
point locations can be computed. Also some of these nodes 
are flagged at this stage so that their displacements can 
be recovered after performing the elastic analysis for 
calculating stress intensity factors by displacement 
correlation method. If only topological information avail- 
able is the connectivity table, this process takes 0(n) 
time where n is the number of elements in the mesh. The 
W-E data structure has a operator (described in 
section 4.3.4) which retrieves a list of vertices adjacent 
to a given vertex in 0(m) time, where m is the number of 
adjacent vertices. 
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4.4.4 Analysis and Remeshing 

Once the geometry and attributes of the problem 
are specified and crack tips identified, the user initiates 
the finite element-boundary element analysis by invoking 
ANALYSIS option in the driver menu page (Fig. 4.14). The 
various important steps in the analysis algorithm for crack 
propagation are as follows: 

1. Perform elastic analysis based upon the geometry of 
the structure, including initial crack geometry, 
appli6<i displacement bomdary conditions and applied 
loads . 

2. Calculate stress-intensity factors and from 
elastic analysis and eqn. (3.44) or eqn. (3.45). 

3. Check wheather a crack should propagate using strain 
energy density criterion (eqn. (2.4)) and calculate 
direction of propagation and amount of propagation, 
also by strain energy density criterion (eqn. (2.11) 
and eqn. (2.4) respectively), if any crack is on the 
verge of propagation. 

4. Update geometry and remesh, if any crack propagates 
by extending structural boundary to reflect the 
propagation . 

5. Return to step 1, With updated geometry and load for 
next step of propagation. 
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The key aspect in modelling discrete crack propagation 
and performing above analysis steps is the ability to modify 
the boundary element-finite element mesh automatically to 
reflect the current configuration. Dozens of man-hour is 
required to produce the results of few crack increments if 
remeshing is done manually [ 71 ] . The program FRAPCO has the 
ability to perform automatic local remeshing around the 
crack tip for each crack increment. The remeshing process 
is performed in following steps: 

1. From the initial direction of crack axis and the 
predicted angle of crack extension with respect 

to crack direction determine the angle of propaga- 
tion in global coordinates. 

2. Replace the quarter-point nodes to their initial 
mid-side position, to remove the local singularity. 

3. Modify the structural boundary to reflect the 
crack extension. 

4. Create a transition mesh that joins the new location 
of crack tip to the tonmodified portion of the 
boundary. 

5. Adjust the mid-side nodes of the new elements around 
crack tip to quarter point location and flag these 
and quarter-points from which stress intensity 
factors will be evaluated in the next step of analysis. 
Display the modified mesh, to allow the user to 
Interactively perform final adjustments. 


6 . 
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An implication of the 3rd step in remeshing, i.e. 
modifying structural boundary, is that when the program 
FRAPCO modifies a certain portion of the mesh, structural 
boundary information must be retained. Thus, before an 
edge is deleted, a check must be made to see if the edge 
makes up a portion of the structural boundary. Since each 
edge is adjacent to two faces, so if one of these faces is 
flagged with the face type indicating that it is not an 
element, the edge is not deleted. The extension of structural 
boundary is then accomplished by the Euler operator SEMV 
(described in section 4.3.3). An existing edge which makes 
up structural boundary is split to extend the boundary to 
describe the new crack configuration. 

The fourth step, i.e. creating a transition mesh to 
join quarter-point elements with the unmodified portion of 
the mesh is achieved by modifying the vertices definitions 
of edges around the crack tip. Edges adjacent to the new 
and old locations of the crack tip are redrawn and the 
modified mesh is displayed for the user to perform minor 
mesh editing before going into next step of analysis. 

In each cycle of analysis, before updating the 
geometry, the user is required to choose between the option 
of updating geometry of all the elements (by adding 
displacements to all nodal point coordinates and local 
remeshing to reflect crack extension) or updating geometiry 
of two boundary element regions around the crack tip (by local 
remeshing to reflect the crack extension only) . In case of 
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material, the finite elements away from the crack 
<ao not deform considerably and thus fairly accurate 
results could be obtained just by updating geometry and 
J^^calculating the two boundary element stiffness matrices 
o^ly. This amounts in substantial saving in cost of analysis 
as finite element stiffness matrices need not be recalculated. 

The stress intensity factors are displayed at the 
end Of second-step of each analysis cycle. User can use any cri 
tsrion . to calculate amount and angle of propagation or 
use the built-in strain energy density function to get the 
same. This type of . interaction of FRAPCO itself can aid 
in the technique and theory development cycle. Because 
the program can at any time display graphically its current 
state and prompt the user for input of certain important para- 
meters (such as angle and amount of propagation) , the fracture 
propagation algorithms need not initially be developed in a 
complete form. Operations which will ultimately be performed 
under algorithmic control can be performed or guided by 
the user manually. This ability not only speeds the technique 
development cycle, it also often adds insight into how 
control algorithms should be constructed, in this way a 
new theory can be tested and accepted or rejected. 

The coupling strategy is employed in the program 
FRAPCO in such a way that it can perform numerical analysis 
by BEM alone or by FEM alone or by a coupled BEM-FEM in a 
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unified >nanner. This kind of integrated fracture mechanics 
analysis system is very useful in the sense that both 
BEM and FEM are available in a single code, so that the 
user can use any one of them at his own discretion. The 
user can start with any kind of mesh and interactively 
change it to put either a few finite elements or a few 
boundary elements around the crack tip at any point during 
the analysis. Thus debate regarding which method is better 
for a particular problem becomes superfluous. 


CHAPTER 5 


RESULTS AND DISCUSSIONS 


5 . 1 Introduction 

In this chapter results of analyses performed with 
the program FRAPCO are presented. The accuracy and conver- 
gence of solution obtained by symmetric bomdary element 
formulation in fracture mechanics problem is studied with 
refining discretization in section 5.2. A comparison of 
time and storage requirements of coupled BEM-FEM and FEM 
is presented next in section 5.3. After getting a feeling 
about convergence characteristics, time and storage require- 
ments the program FRAPCO is tested with a few test problems 
representing different kinds of fracture mechanics problems 
in section 5.4. Having checked the working and accuracy of 
the program, a stable crack propagation study of a centre 
crack in a rectangular plate is undertaken in section 5.5. 
The effect of initial crack size on various parameters 
such as propagation initiation load, critical load for 
bringing the crack on verge of global instability etc. is 
also presented in section 5.5. Finally, section 5.6 
] 3 j-;j_ngs out the results of mixed mode crack propagation of 
an angled notch. 
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^ • 2 Convergence of Symmetric BEM ForTnulation 

The success of a coupled BEM-FEM analysis depends 
on the accuracy of symmetric boundary element formulation. 

In order to assess the accuracy obtainable by the symm- 
etric coupling theory, a numerical convergence study is 
undertaken with the aid of a mixed— mode centre slant crack- 
problems, as shown in Fig. 5.1. The accuracy of stress 
intensity factors and are studied since their correct 
determination ensures correctness of calculation of the 
amount of propagation and angle of propagation determined 
through equation (2.4) and equation (2.11) respectively and 
also because they reflect the correctness of the stress- 
analysis near the crack tip. 

The problem is analysed with fixed size singular 
elements (1/a =0.10 where 1 = length of crack tip boundary 
element and a = half crack length) at the two crack tips 
and increasing number of non-singular quadratic boundary 
elements over the rest of domain. The problem is represen- 
tative of a crack in an infinite domain. The results are 
compared with the exact results of a crack in an infinite 
plane [75] : 

t — 2 

Kj =o/Tra Sin $ 

Kjj = o/ir^ Cosg Sing 

» 

in which g is the angle between the vertical load axis 



cr - 1 ksi 


t t t t t t t f f f f f f 



E = 5250 ksi , 1>=0'2 . W = 10", 
2a = >/F inch , p = 45 



Crack lip element 

FIG. 5*1 Typical BEM hfiesh for centre slant crack problem 



(excluding crock tip elements) 

FIG. 5-2 Percent error in the computed SIF for slont crock problem 
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and the crack axis as shown in Fig. 5.1 . No finite 
length and width correction factors are applied to the 
SIFs of infinite plane problem as these correction 
factors would be quite small [75]. 

Results of the convergence study is presented in 
Fig. 5.2, in which results of symmetric formulation are 
plotted along with the results of conventional BEM [18 ] 
with same number of subdomain elements. From Fig. 5.2, 
it is clear that a better rate of convergence is obtained 
by conventional BEM as indicated by the steeper slopes of 
curves. But even though slower rate of convergence is 
achieved by the symmetric BEM formulation, a reasonable 
number of subdomain elements produce results comparable to 
that obtained by conventional BEM. 

5. 3 A Comparison of Time and Storage Requirements of 
Coupled BEM-FEM and FEM 

For comparing the coupled BEM-FEM with FEM from 
time and storage requirements view point, a centre crack 
problem as shown in Fig. 5.3 has been chosen. The 
discretization for the coupled method is obtained by 
removing few finite elements from the crack tip creating 
a botondary element region there. For comparison, the 
number of nodes on boundary is kept same for both 
FEM and coupled method. Also, in both the cases, singular 
quarter point elements are used around the crack tip, as 



F16. 5-3 Centre crack problem (a) Finite element mesh 
(b) Coupled 8EM-FEM mesh 



shown in Fig. 5.3. Similar interactive graphics prepro- 
cessors are used to input both the meshes and all compu- 
tations are done on same ND-560/CXA machine. 

Table 5.1 summarizes the salient characteristics 
of time requirements for FEM and coupled method. As 
observed from table 5.1, coupled method performs better 
than FEM in terms of CPU times with comparable accuracy 
of results being obtained. The equation formulation 
time is less because less number of nodes and elements 
are used in coupled analysis for which less computation 
and less communications between secondary disc storage 
units and core memory are required. Equation solution 
time is also less because of resulting smaller system of 
equations and also because the Frontal solver has to 
address the secondary disc storage units less number of 
times. For the coupled method, the preprocessing and 
other related operations, called data-handling time in 
Table 5.1, shows marked improvement over FEM in time require- 
ments. This is because less number of edges, faces and 
vertices are used in coupled boundary element-finite element 
mesh resulting in a substantial reduction in mesh display, 
mesh editing, attribute editing and remeshing timings. As 
pointed out in section 4.3.5, the time complexity of 
topological queries (TOl of W-E data structure is 0(6EV^+3K) 
where K is constant access time and EV^ denotes the linear 



TABLE 5.1 : Centre Crack Problem: Comparison of 

Typical PEM and Coupled BEM-FEM 
Solution Times 



PEM 

Coupled 

BEM-FEM 

Preprocessing time 
(man-minute) 

6 

6 

Number of Elements 

40 

9^ 

Number of Nodes 

149 

76 

Equation formulation time 
(CPU-sec)* 

243 

132 

Equation solution time 

135 

60 

(CPU-sec) * 



Data hand] ing time 

* 

(CPU-sec) 

51 

29 

Total solution time 
(CPU-sec) * 

1 

1 

469 

221 

Mode I stress intensity 
** 

factor 

1.190 

1.189 

(*) nD-IOO+ND- 500 CPU-Time 

9 



(**) Ref. [83]- 1.187. 

(#) Number of boundary elements in boundary 
element region - 24 



















118 


dependence of access time on the local quantity EV^, 
the number of edges around vertex V^. As the quantity 
EV^ is much less in coupled boundary element-finite element 
mesh than finite element mesh, so in general, access 
time for topological queries are less in coupled method. 

As pointed out in section 4.4, the TOs are used hundreds 
of time during mesh generation, mesh editing, remeshing, 
attribute editing, preparation of data compatible to 
numerical analysis subroutines (e.g. the VF^ query is 
used for each element to extract the node-element connecti- 
vity from W-E data base during stiffness matrix calculation) . 
This advantage would be more prominent in analysis of 
large system where thousands of edges, faces and nodes 
are used, making the data handling time comparable to the 
numerical solution time. 

In this connection it is interesting to note the 
time and algorithmic complexity involved in the remeshing 
to reflect configuration change due to crack propagation. 

In program FRAPCO boundary elements are used around the 
crack tip which makes the algorithm of remeshing (described 
in section 4.4) considerably simpler than remeshing 
algorithms of FEM, e.g. as described in [47 ] . As only 
three quarter point nodes are used around the crack tip 
node, so. the search made by query to locate such 
nodes around the crack tip takes much less time than 
corresponding time for FEM remeshing. In [47 ] , a remeshing 
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Strategy is described in which a number of VF. and EF. 

1 1 

queries are made which is not required in the remeshing 
algorithm used by FRAPCO, making amomt of overall 
remeshing time for each crack much smaller. 

A comparison of the storage requirement for the full 
topological data structure of the meshes of Fig. 5.3 is 
shown in Table 5.2. Assumptions made to allow comparison 
are that all the fields of each record structure are 32 bits 
in length and that an eight-bit byte is the minimum-size 
storage unit. Sizes of face (40) , vertex (24) and edge (40) 
records were derived from the support record structure defined 
in section 4.3.8, All the fields of each record structure 
were considered, though some of the fields are implemen- 
tation dependent, e.g., the maximum and minimum values of 
X- and Y-coordinate values of a face could be dropped in 
a storage and time trade-off. These choices have a tendency 
to minimize the apparent difference between FEM and coupled 
BEM-FEM storage requirements; the figures presented are 
therefore sort of best-case difference . An implementation 
which tries to reduce storage requirements by dropping some 
of the fields (which could be accessed through some search 
algorithms thereby increasing the time complexity) would 
increase the percentage of space attributed to differences 
in the FEM and coupled BEM-FEM meshes. 









121 


For the type of imesh considered, finite elements 
mesh requires about 56% more storage than coupled BEM-FEM 
mesh. This is to be expected as finite element mesh 
consist of more number of edges, faces and vertices. As 
pointed out in section 4.3.5, considering only the pointers 
absolutely necessary, the W-E data structure storage require- 
ment is 9E where E is the total number of edges in the 
graph. This means that storage requirements increases 
linearly with E. As FEM mesh has a higher edge to face 
ratio than the coupled BEM-FEM mesh, the storage requirement 
for FEM shows dramatic increase. This factor is of utmost 
importance, particularly in case of very large structure 
where thousands of edges, faces, vertices are used and it 
tilts the balance in favour of coupled analysis. 


5.4 Testing of the Program 

Having obtained a fair idea about the number of 
elements to be used to get acceptable results from the 
convergence study presented in section 5.2 and associated 
storage and time complexity requirements from section 5.3, 
testing of the program FRAPCO is undertaken in this section. 
Three fracture mechanics problems are studied in which 
stress intensity factors of a given cracked body are 
determined and compared with analytical and/or numerical 
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results reported in literature. Again the stress intensity 
factors are taken for the purpose of comparison as they give 
fair estimates of accuracy of the solution obtained by FRAPCO, 
as pointed out in section 5.1. The material values used for 
the problems in this section are E = 5250 ksf and v = 0.20. 


TEST PROBLEM 1 

The coupled BEM-FEM is applied to solve a centre 
crack problem. A typical mesh shown in Fig. 5.4. The 
problem is doubly symmetric about the crack axis and hence 
only one quarter of the plate is considered for analysis. 
The mode -I stress-intensity factor is computed for a fixed 
height to width ratio (V-b) of 2.0 but varying half crack 
length to width ratio (aA») • The results are presented 
in Table 5.3 and they are compared with that of Isida [83 ]. 

TEST PROBLEM 2 

The branch crack problem of Fig. 5.5 is studied 
next. This kind of branching occurs during many situations 
in real life such as, mixed mode loading, dynamic crack 
propagation etc. The problem, shown in Fig. 5.5, is 
composed of two branches of equal length in a square plate 
under uniform tension. A very large value of W/a (= 50) 
is taken and in Table 5.4, stress-intensity factors at tip 
A are compared with that of Kitagawa et.al. [ 80 3 who solved 
the same problem for an infinite plate. 



FIG. 5-5 Branched crack problem 
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TABLE 5.3 : Stress Intensity Factor for a Centre 

Crack: Test Problem-1 


a/b 

K 0 /trs 

Kj/cj /ira* 


-present 

-Ref [ 83 3 


analysis 


0.3 

1.068 

1.058 

0.4 

1.126 

1.109 

0.5 

1.208 

1.187 

0.6 

1.314 

1.303 


TABLE 5.4 : Stress Intensity Factor for a Bent 
Crack : Test Problem-2 


2ira(l+CosBl 


Kj /2ira(l+Cosej 


Present 
analysis 
(W/a = 50) 


Ref. [80 ] 
(W/a = ») 


Present 

analysis 

(w/a=50) 


Ref [so 3 

(W/a = «) 


0.574 


0.596 


0.654 


0.641 
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TABLE 5.5 ; Stress Intensity Factors for a Crack 

Emanating from a Circular Hole: 

Test Problem- 3 


a/R 


0.2 

0.3 

0.4 


K^/af~2V 

Kj/a /2T 

Present 

Ref. [ 81 ] 

analysis 

(W/a=50) 

(W/a=*) 

2.36 

2.30 

2.11 

2.04 

1.93 

1.86 


0.5 


1.85 


1.73 
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TEST PROBLEM 3 

A crack emanating from a circular hole in a plate 
subjected to uniaxial tension (Fig. 5.6) is solved next 
by the present program FRAPCO. The problem represents 
a crack in a region of high stress concentration. A 
typical mesh to perform the analyses is shown in Fig. 5.6. 
Again a very large value of W/a (=50) is taken and results 
are compared with Bowie's results [si,}, who solved the 
problem for an infinite plate. Table 5.5 brings out the 
comparison of mode I stress intensity factor values of crack 
tip for various hole radius. 


5.5 Propagation of a Centre Crack 

The process of slow stable growth of a centre 
crack in a rectangular sheet is studied using the present 
program FRAPCO. A plate of height 2h and width 2b with a 
centre crack of length 2a (Fig. 5.7) is considered. The 
plate is subjected to uniform \iniaxial tensile stress 
perpendicular to crack plane. The influence of crack's 
geometry on the process of stable crack growth is also 
addressed and results are compared with that of Gdouto [79], 
The width 2b and height 2h of the panel are taken to be 
20 inch and 16 inch respectively, while the crack length 
is varied. 





FIG. 5-7 Center crack problem 
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The crack propagation initiation load, is 
calculated using strain-energy density criterion for 
various initial crack lengths. Fig. 5.8 presents vari- 
ation of versus the initial crack length. As seen 
from the Fig. 5.8, decreases with initial crack length. 
Thus longer cracks require less initiation load as stress 
intensity factors increase rapidly with the crack length. 

After the crack initiation, the crack grows in a 
stable manner until global instability is reached. In 
order to study the process of stable crack growth, the 
applied stress a is increased by constant interval 
and the corresponding crack increments Ar are calculated 
by strain-energy density criterion, as explained in section 
2.3 The process is continued until the crack extension 
size Ar becomes equal to critical ligament size r^, which 
corresponds to global instability. Table 5.6 presents the 
successive crack increments Ar of a growing crack of 
length a = 2 inch, as stress a is increased at a constant 
interval of Aa = 2.0 ksi . The first value of the applied 
stress in the table = 38. 7 ksi corresponds to the 
critical stress for crack initiation, and the last value 
= 64. 7 ksi corresponds to the critical stress for 
global instability. The corresponding results of Ref. [79] 
are : a. = 39. 8 ksi and a = 58. 2 ksi. Though crack 

X c 

propagation initiation load corresponds well with that 
of Ref. [79] , the critical stress corresponding to global 
instability shows significant deviation from that in 
Ref. [ 79 ]. This is expected because in Ref. [79] , load 


TABLE 5.6 


Values of Relevant Parameters for Various 
Crack Steps of a Growing Centre Crack 


Number of 
crack incre- 
ment 


Applied 
stress 
a (KSi) 


Crack 

increment 

r(in)xlO”^ 



38.7 


40.7 


42.7 

44.7 

46.7 

48.7 


50.7 


52.7 


54.7 


56.7 


58.7 


60.7 


62.7 


64.7 


1.0816 

1.1937 

1.3116 

1.4351 

1.5644 

1.6996 


1.8405 


1.9874 


2.1403 


2.2991 

2.4641 


2.6352 


2.8126 


2.9963 



2.001081 

2.002275 

2.003586 

2.005021 

2.006585 

2.008285 

2.010125 

2.012113 

2.014253 

2.016552 

2.019017 

2.021652 

2.024465 

2.027461 








FIG. 5-8 Critical stress for crack initiation vs. initial half 
crack length 



FIG. 59 


Critical stress for global instability vs. Initial half 
crack length 





FIG. 5-10 Variation of difference of critical stress for insta 
bility and crack initiatian ,<rc - (Tj ,ys. Initial 
half crack length 
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steps taken is La — 0.2 ksi and crack was propagated 
by 93 steps to arrive at the Instability point whereas 
in the present analysis lead steps Aa= 2.0 ksi and 
number of propagation steps 14 are taken. This shows 
the strong load step dependence of stable crack growth. 

Figure 5.9 shows variation of stress a versus 

c 

the initial crack length a. It is observed that o 

c 

decreases as a Increases. In order to obtain an estimate 

of the stress increment from crack initiation to unstable 

crack extension, variation of (a - a ) versus the initial 

c i 

crack length a is presented in Fig. 5.10. It is observed 
that stress increment needed to bring the crack tip to 
onset of global instability from initiation of propagation 
decreases with initial crack length. 


5.6 Propagation of an Angled Notch 

The quasi-static crack propagation analysis of an 
angled notch problem shown in Fig. 5.11 is studied in 
this section. The present analysis results are compared 
with the same analysis using finite element method as 
well as experimental results of Ingraffea [7?] . The 
initial mesh used and the bovindary conditions considered 
for coupled BEM-FEM analysis is shown in Fig. 5.13. 

The material is Indiana (Salem) limestone the properties 
of which are indicated in Fig. 5. 11. 



w 



W r 4‘*,t = 3/4". 
a s.4",b =.004*: 

E = 5*25 X lO^psI, 

V = 0.21,ro=.00l5'; 
(dw/dv)cr =0-39in-lb/in3, 
K|c= SOOPslVTn 


FIG. 5 11 Angle notch plate, edge loaded in 
uniaxial compression 



FIG. 512 Location of fracture Initiation point 
on surface of ellipse 
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Since the problem is that of a notch, the surface 
layer energy concept [ 50] is used to predict the location 
for crack initiation point. The notch is assumed to behave 
like an idealized ellipse f 78] , with an effective semi -minor 
axis computed from. 




(5.1) 


where a = notch half-length, 0.400 inch, and P= notch tip- 
radius, 0.004 inch. The point of maximum tensile tangen- 
tial stress p^ is computed according to the positive root 
solution of 


tan 


b [a Sin^ g 


a 


-b Cos^g + Va^ Sin^g+ b^Cos^g ] 
(a+b) SinB Cos 6 


(5.2) 

The geometric relationship between the eccentric angle 
and p^ is depicted in Fig. 5.12. The values of x-coordinate 
of p^, and r\^ for different notch angles are listed in 

Table 5.7. 


Fracture initiation load, p^, is calculated according 
to surface layer strain energy density theory for notch [ 50 ] 
which states that fracture initiation occurs when (‘SW/dV)^^^^^ 
reaches a material critical value, (dW/dV)^^ at a distance r^ 
from the point Pj, tx^, ,ti£)[ as shown in Fig. 5JL2], the fracture 
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initiation point on the notch surface . The initiation 
load versus notch angle g is then obtained by norinalising 
this load according to, 

^i f ^ = F (B) (5.3) 

The results obtained by present analysis are compared 
with experimental values obtained by Ingraffea [82 ] in 
Fig. 5.14. As observed, the results agree quite well with 
the experimental data for $ ^ 45°. This is expected because 
for B > 45°, notch is a poorer approximation to the ellipse 

for predicted in this region. Also, as pointed out by 
Ingraffea [82], the area over which the tensile stress concen- 
tration acts increases rapidly with 3 . This implies that a 
slight local material weakness, such as grain bomdary or weak 
grain, located any where within such a region may trigger 
fracture initiation well away from the theoretical point, 
resulting a wrong prediction of initiation load value. 

The program Frapco is next used to perform crack 
trajectory prediction analysis. The notch angle B is taken 
to be 45°. The initial mesh used to perform the coupled 
BEM-FEM analysis is shown in Fig. 5.13. The load is 
gradually increased in steps and crack path is predicted 
from strain energy density theory, as outlined in Chapter 2. 
The computed crack path is compared with the experimentally 
obtained results and FEM results [77] , in Fig. 5.15. The 





1.77 

2 


o FRAPCO 

A FEM [77] 

• EXPTL,LOWER 
NOTCH [77] 

P in ksl 



FIG. 515 Comparison of present analysis crack 

path with exptl, results and FE analysis 
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FEM results are obtained by maximum strain energy release 
rate theory. As seen from Fig. 5.15, there is good agree- 
ment between present analysis and experimental results. 

However, the finite element crack path prediction, parti- 
cularly the load at certain steps, agree less with the 
present analysis. This discrepancy is most likely due to 
small errors in crack topology introduced by the irregularity 
of a finite element mesh. Using boundary elements aroxind crack 
tip, the crack trajectory is modelled more smoothly in the 
present analysis and element distribution is also more regular. 
The resulting stress factor evaluation is thus more accurate, 
which in turn, results in a better load prediction at each 


crack step. 



CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 

6, 1 Conclusions 

The results of the present study have shown that an 
interactive discrete crack propagation analysis code based on 
a coupled BEM-FEM is a very efficient technique. The follow- 
ings are the specific conclusions of the present work; 

1. The coupling of BEM and FEM exploits the advantages of 

both methods. By using boundary elements aroTjnd the 
crack tip, the following advantages are gained - 
(a) better accuracy of solution, (b) ease of discre- 
tizing a moving boundary, (c) ease of the use of 
singular elements, (d) considerable saving in computer 
time due to less number of nodes and elements as 
boundary needs to be discretized. On the other hand, 
use of FEM away. . from crack tips offer the following 
advantages - (a) better efficiency of FEM for analysis 
of finite bodies with large dimensions, (b) use of 
locally based interpolation fiinctions which produce 
narrowly banded final set of equations, (c) avail- 
ability of frontal solution technique which requires 
less core storage and has got resolution facility. 



Considerable saving in both time and storage require- 
ments can be achieved by using a coupled BEM-FEM 
compared to FEM. This, combined with the fact that a 
unified approach can be adopted in finite element 
and boundary element regions mesh generation, data 
base design, stiffness matrices’ assembly etc, make 
such coupled BEM-FEM an attractive choice as it could be 
incorporated . into an existing interactive finite 
analysis system with minor modifications. 

Ability of a program to perform numerical analysis by 
BEM alone or FEM alone or by a coupled BEM-FEM is 
very advantageous as the user can use any one of these 
methods at his own discretion. 

Transfinite mapping mesh generation technique could be 
used to generate mesh for finite element region and 
bomdary element region equally efficiently. 

The winged edge data structure is very efficient for 
interactive finite element and coupled boundary element - 
finite element analyses because - (a) All topological 
adjacency queries can be satisfied in constant or linear 
time, (b) data base can be modified to reflect topology 
change as a result of crack propagation without 
reorganizing mmodified portion of the data, (c) program 
response time grows at a rate no more than linearly 
with the problem size, (dl higher level operator could 
be defined which hide the Inherent complexity of creating 
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and modifying this type of data structure. 

6. Integrating all the aspects of a fracture mechanics 

analysis system such as fracture propagation theory, 
boundary and finite element methods, mesh generation, 
data base design etc. in a single code can aid in the 
technique and theory development cycle. Any new 
fracture propagation criterion and algorithm could be tested 
and accepted (or rejected) very efficiently and easily. 


6 . 2 Recommendation for Further Study 

The computer code developed herein, though quite versatile 
and general, many improvements over the present algorithms can 
be made. Also modular nature of the program leaves the possi- 
bility of incorporating new algorithms wide open. Few suggestions 
towards the development of a comprehensive fracture mechanics 
code by working on the present program are enlisted in the 
following lines. 

1) The program FRAPCO relies heavily on virtual memory 

capability of the machine. A data structure, such 
as in [60 3 to take into account the restriction 
on the in-core storage, easy checking and access to 
each of the arrays involved and use of auxiliary memory 
could be employed to rm the program in a medium size 
non-virtual memory machine. 
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2) In corporation of adaptive numerical integration 
scheme [ 60 ] in boundary element formulation is 
suggested. 

3) Special crack tip elements [66] , [84 jother than 
quarter point elements could be used to improve 
the accuracy of the results. Particularly attra- 
ctive is the Hermitian cubic singular elements [ 66 ] 
as they do not require any division of problem 
domain into subdomains to model the crack surface. 
Thus multidomain BEM formulation, which requires 
additional effort from the user and also introduces 
error in the solution by constricting the solution 
along the interface to conform to the functional 
variation of the elements, can be avoided. 

4) Extension to elastic-plastic analysis with few finite 
elements around crack tips and boundary elements 
away from it could be achieved very easily. 
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